## 3 Methods
### Spatial-Temporal Extreme High XCO\({}_{2}\) Detection  
We define extreme high XCO\({}_{2}\) as the spatiotemporal continuum units whose fitting residuals deviate outside of the normal range of high probability values [7; 30]. For each grid, we calculate a Z score of fitting residuals to standardize across differences in the magnitude of residuals according to Equation (2). The Z score, \(Z(\mathrm{i},\mathrm{j},\mathrm{t})\), is calculated from the XCO\({}_{2}\) fitting residuals, \(\delta\) but also described as dXCO\({}_{2}\) throughout the manuscript, in a grid with latitude i, longitude j, and time t, and mean\((\delta\;(\mathrm{i},\mathrm{j}))\) and std\((\delta\;(\mathrm{i},\mathrm{j}))\) represent the mean and standard deviation value of the time series in location \((\mathrm{i},\mathrm{j})\).  
\[Z(\mathrm{i},\mathrm{j},\mathrm{t})=(\delta(\mathrm{i},\mathrm{j},\mathrm{t})- \mathrm{mean}(\delta(\mathrm{i},\mathrm{j})))/\mathrm{std}(\delta\;(\mathrm{ i},\mathrm{j})), \tag{2}\]  
We selected the Gaussian model, a simple spatiotemporal outlier detection technique to define anomalously high XCO\({}_{2}\) in Equation (3a). This test identifies grids where dXCO\({}_{2}\) (or \(\delta\)) is much larger than a threshold, \(\mu_{1}\), and has a Z score of the dXCO\({}_{2}\) variability that exceeds a threshold, \(\mu_{2}\). In this study, \(\mu_{1}\) is set to be 1.0, according to a priori knowledge of XCO\({}_{2}\) concentration's precision from GOSAT retrievals, and \(\mu_{2}\) is set to 1.96, which represents the 95% threshold in the Gaussian distribution. The logic test is described using 'bool'; true = 1 and false = 0. The connectivity of the extreme units are tested as well. D\((\Delta\mathrm{i},\Delta\mathrm{j},\Delta\mathrm{t})\) represents the spatiotemporal distance between grids (Equation (4)). The distance threshold is \(d\) and the connectivity threshold is n. Thresholds d and n were set to be 1 (minimum distance) and 3 (mid-value of range from 0 to 6), respectively. In order for an event to count as an extreme event, Equations (3a)-(3c) must all be true.  
\[\mathrm{bool}\big{(}\delta\;(\mathrm{i}_{0},\mathrm{j}_{0},\mathrm{t}_{0})> \mu_{1}\big{)}\times\mathrm{bool}\big{(}Z(\mathrm{i}_{0},\mathrm{j}_{0}, \mathrm{t}_{0})>\mu_{2}\big{)}=1, \tag{3a}\] \[\sum_{\mathrm{i}=\mathrm{i}_{0}-\mathrm{d}}^{\mathrm{i}=\mathrm{ i}_{0}+\mathrm{d}}\sum_{\mathrm{j}=\mathrm{j}_{0}-\mathrm{d}}^{\mathrm{i}= \mathrm{j}_{0}+\mathrm{d}}\sum_{\mathrm{t}=\mathrm{t}_{0}-\mathrm{d}}^{\mathrm{ t}=\mathrm{t}_{0}+\mathrm{d}}\{\mathrm{bool}(\delta(\mathrm{i},\mathrm{j},\mathrm{t})>\mu_{1}) \times\mathrm{bool}(Z(\mathrm{i},\mathrm{j},\mathrm{t})>\mu_{2})>\mathrm{n}\},\] (3b) \[\mathrm{Where}\;\mathrm{D}\big{(}\mathrm{i}-\mathrm{i}_{0},\mathrm{j}- \mathrm{j}_{0},\mathrm{t}-\mathrm{t}_{0}\big{)}<\mathrm{d},\] (3c) \[\mathrm{D}(\Delta\mathrm{i},\Delta\mathrm{j},\Delta\mathrm{t})=\sqrt{( \Delta\mathrm{i})^{2}+(\Delta\mathrm{j})^{2}+(\Delta\mathrm{t})^{2}}, \tag{4}\]  
In addition, three-dimensional extreme dXCO\({}_{2}\) can be inspected with the steps above. We defined the extreme XCO\({}_{2}\) events as the space-time continuous unit, e.g., if one grid is selected as extreme high dXCO\({}_{2}\), there are six directions for expansion. If any grid, within the six directions, is also an extreme high dXCO\({}_{2}\), they are included in the same unit. They can further expand in time as well, until no other grid are connected to the unit and we count the total number of grid units to determine the magnitude of the extreme event. Identified units of the space-time continuous grids are then the extreme high XCO\({}_{2}\) units. This procedure is the flood-fill algorithm described in [PERSON] et al. (2017). We used the function of bwconncomp from Matlab to execute the algorithm described above. The work flow of the extreme events detection method is summarized in Figure 1.