**1. Introduction**

Global Navigation Satellite System Reflectometry (GNSS-R) data have the potential to regionalize methane (CH4) emissions from land surface images by detecting their inundation status. Methane is an important greenhouse gas (GHG); its global warming potential over a 100-year horizon is 28 times higher than that of carbon dioxide (CO2) [1]. In 2011, the CH4 concentration was 1803 ppb, 150% higher than the preindustrial level, and a predominantly biogenic post-2006 increase has also been reported [2]. Concurrently, atmospheric methane's *δ*13CCH4 value has trended towards lighter (13C-depleted) values, implying a significant shift in the balance between the sources and sinks of CH4 [3] and a greater contribution of biogenic CH4 emission sources rather than fuel combustion to this rapid CH4 concentration increase [2]. Several hypotheses have been postulated for the cause of this isotopic shift, and these hypotheses can be summarized as one or a combination of the following: (i) a change in the oxidative capacity of the atmosphere [4]; (ii) changes in the relative strengths of anthropogenic sources, such as land-use changes on tropical wetlands to agriculture or waste and fossil fuel emissions with an overall net effect

of increasing emissions (e.g., [5]); and (iii) an increase in natural sources such as wetlands, potentially as a feedback effect from regional climatic change (e.g., [3]). Large gaps still exist between top-down and bottom-up CH4 total global emissions calculations, with much of the uncertainty associated with the emissions of wetlands and other natural emissions categories [6,7], particularly in tropical wetlands [7–9].

Because CH4 is emitted from inundated soil, which is spatiotemporally heterogeneous and has a flux pattern characterized by non-Gaussian/nonlinear behaviors [7], the appropriate evaluation of the CH4 flux requires the monitoring of the inundation status with spatiotemporally high-resolution techniques [10]. GNSS-R data became a popular input source in microwave remote sensing techniques following the deployment of the Cyclone Global Navigation Satellite System (CyGNSS), an eight-microsatellite constellation data system [11]. Every single CyGNSS microsatellite has two left-hand circular polarization (LHCP) down-looking antennas pointing to the Earth's surface with an inclination angle of approximately 28 degrees on either side of the satellite ground track [12].

The data can be used to globally detect the land surface inundation status almost daily with high-spatial-resolution L-band microwave signals (with estimated spatial resolutions of approximately 500–7000 m [13]) compared to common passive L-band microwave radiometers. A few studies have reported that the use of CyGNSS-based inundation maps for land surface methane emission simulations improved the representation of the CH4 emission status compared to the results obtained using common wetland maps (e.g., simulating a greater amount of CH4 emissions by detecting inundation under clouds/vegetation with GNSS-R data [14]).

There are several studies on the detection of inundation over wetlands with GNSS-R data e.g., [15–19]. However, the results in most studies remain spatiotemporally sparse. In most cases, the spatiotemporal interpolation is conducted with monthly observation datasets, or spatially interpolated with optical observation sensors e.g., [15,19]. Due to the limitations of L-band fine-spatial-resolution microwave remote sensing data like GNSS-R, there are only a few studies conducting the cross-validation of GNSS-R and L-band SARs observations [15]. Furthermore, from the perspective of the application of this study, most of the time, this sort of fine-spatial-resolution, satellite-derived wetland/inundation observation is downsampled or spatially thinned (a.k.a., superobservations) before being used in advanced simulation modeling approaches accompanied with high computation costs (e.g., coarse-spatial-resolution ensemble simulations or the use of superobservations to deal with observation error covariance in data assimilation tasks) by degrading the spatial resolution or thinning the observations (e.g., [14,20]). Due to the local heterogeneity of the inundation status and the non-Gaussian/nonlinear characteristics of the spatiotemporal CH4 emission distribution at the local scale [7–9], the deterioration of the spatial resolution of data can introduce large discrepancies to the emission values obtained between the top-down approach and bottom-up approach [6–9,20]. Therefore, the regionalization of CH4 emissions based on high-spatial-resolution L-band microwave data as a bottom-up approach still remains important [7–9,17]. Since most studies have used GNSS-R data for regional-scale simulations at a relatively coarse spatial resolution compared to remote sensing observations (e.g., 0.01◦-resolution CyGNSS-based watermasks are downsampled to a 0.5◦ resolution to match the WetCHARTs simulation grid [14]), few studies have paid attention to the differences among each specular point's footprint size (i.e., the glistening area). To rasterize each piece of specular point-scale vector data without downsampling for use in local-scale simulations, one must consider the difference among each specular point's footprint size to use these signals in fine-spatial-resolution, local-scale simulations (e.g., 10–50 m resolution irrigation models [10]). This information would also be essential for determining the spatial localization scale to ensure efficient data assimilation by determining the spatial localization scale at each specular point and adequately addressing the spatial observation error covariance.

More fundamentally, the amount of data of a certain quality provided by the GNSS-R microsatellite constellation is still limited, and the observations are prone to being conducted sparsely in space; in addition, the incidence angle varies widely among specular points, which is known to cause biases in the microwave reflectivity observations (unlike other spatially continuous microwave remote sensing observations, such as those obtained from passive microwave radiometers or synthetic aperture radars). The local incidence angles of the Phased-Array L-band Synthetic Aperture Radar-2 (PALSAR-2) ScanSAR instruments vary from 25 to 50 degrees, while CyGNSS incidence angles vary from 0 to 70 degrees [13,17]). To prepare inundation maps based on GNSS-R data for future applications or to be assimilated into simulation models, the spatiotemporal interpolation step needs to be processed before the data can be used in applications. Therefore, an adaptive quality control method that considers the size of each specular point and depends on each specular point vector and incidence angle but does not require ad hoc parameter tuning or region-specific empirical parameterization with external data, such as the normalized differential vegetation index (NDVI) or digital elevation model (DEM) data, is essential for this robust interpolation preprocessing step. To realize this globally consistent rasterization at a fine spatial resolution, the authors have developed a precision index calibration scheme implemented while processing the raw specular vector data to rasterize data while considering the differences in incidence angles and specular points' sizes/shapes/velocities. To compensate for the spatially sparse distribution of GNSS-R specular data, the temporal Kalman smoother is applied by using the precision index as the reciprocal observation error number in each 15-day cycle over the Mekong Delta as a demonstration; this case study area consists of double-/triple-rice-cropping systems, aquacultural ponds, mangroves and peatlands. Cross-validation with the PALSAR-2 quadruple polarimetric data (3–6 m resolution) product is also conducted, and the results are validated with ground inundation observation datasets [7–9,17]. The goal of this study is to demonstrate the usefulness of this quality control method by applying it to fine-spatiotemporal-resolution analyses over the Mekong Delta [i.e., (I) comparing it with a common change detection algorithm with the daily temporal resolution, (II) applying it with a 500 m rasterization with a 15-day temporal resolution, (III) and surveying the consistency with 3–6 m spatial-resolution L-band SAR backscatter intensities].
