Gap‐Filled Multivariate Observations of Global Land–Climate Interactions

The volume of Earth system observations has grown massively in recent decades. However, multivariate or multisource analyses at the interface of atmosphere and land are still hampered by the sparsity of ground measurements and the abundance of missing values in satellite observations. This can hinder robust multivariate analysis and introduce biases in trends. Nevertheless, gap‐filling is often done univariately, which can obscure physical dependencies. Here, we apply the new multivariate gap‐filling framework CLIMate data gapFILL (CLIMFILL). CLIMFILL combines state‐of‐the‐art spatial interpolation with an iterative approach accounting for dependencies across multiple incomplete variables. CLIMFILL is applied to a set of remotely sensed and in situ observations over land that are central to observing land–atmosphere interactions and extreme events. The resulting gridded monthly time series covers 1995–2020 globally with gap‐free maps of nine variables: surface layer soil moisture from European Space Agency (ESA)‐Climate Change Initiative (CCI), land surface temperature and diurnal temperature range from Moderate‐resolution Imaging Spectroradiometer, precipitation from GPM, terrestrial water storage from GRACE, ESA‐CCI burned area, and snow cover fraction as well as 2‐m temperature and precipitation from CRU. Time series of anomalies are reconstructed better compared to state‐of‐the‐art interpolation. The gap‐filled data set shows high correlations with ERA5‐Land, and soil moisture estimates compare favorably to in situ observations from the International Soil Moisture Network. Soil moisture drying trends in ESA‐CCI only agree with the reanalysis product ERA5‐Land trends after gap‐filling. We furthermore showcase that key features of droughts and heatwaves in major fire seasons are well represented. The data set can serve as a step toward the fusion of multivariate multisource observations.

States of America are aiming to observe essential climate variables from space, wherever possible.However, the most prevalent approach for these activities is to focus on each variable separately, which can obstruct further analysis.For this reason, the need to fuse these data into coherent representations-also referred to as "Digital Twins"-of the Earth has been voiced by many actors in the field of Earth observations (Bauer et al., 2021;Mahecha et al., 2020;Salcedo-Sanz et al., 2020).
Moreover, assessing cross-variable dependencies and tracking multivariate (or compound) extreme events with a multisource or multiobservation analysis is still hindered by a high fraction of missing values in many remotely sensed or ground observations.Satellite observations often have ample missing values due to retrieval limitations or flight geometry (Bessenbacher et al., 2022).Ground observations are sparse and have insufficient coverage in many regions in the world (Bessenbacher et al., 2023).Therefore, combining two or more incomplete observations remains difficult, as the different patterns of missing values (also referred to as "missingness") often do not overlap; therefore, only a small fraction of points can be retrieved where both variables are observed.These large-scale and ubiquitous gaps in observations are unavoidable: for example, in soil moisture products produced by the ESA-CCI, only around 40% global land data points are observed in the years 2003-2020, which can be an obstacle for several applications including large-scale drought monitoring or the investigation of processes at the land-atmosphere interface (Bessenbacher et al., 2022).
For this reason, implicit or explicit gap-filling is common.First, many freely available data products have been subject to implicit gap-filling.For instance, the Moderate-resolution Imaging Spectroradiometer (MODIS) land surface temperature (LST) data set ignores gaps in daily observations, that is, the monthly averages are only based on the observed days, which might be different across regions (Wan et al., 2015b).ESA-CCI burned area (BA) is available in monthly resolution, with the fraction of observed area per grid point per month given, which does not allow for recovery of the original missing value pattern of the satellite retrieval (Otón, 2020).Similarly, bottom of atmosphere fluxes in the ESA-CCI radiation data set count the number of observations upon which a monthly value is created, but the times of these observations are not available (Hollmann & Pinnock, 2020;Stengel et al., 2020).Second, many studies employ explicit gap-filling prior analysis to use incomplete observations from multiple sources.These treatments usually build upon spatial or temporal interpolation and are on a spectrum from simple gap-fill with a mean seasonal cycle (Tramontana et al., 2016), spatiotemporal interpolation (Zhang et al., 2022) to sophisticated Kalman-smoother time interpolation (Huffmann et al., 2019).Third, explicit gap-filling frameworks that exploit spatial, temporal, or spatiotemporal relations are common for popular remotely sensed observations from the MODIS instrument (Bai et al., 2023;Boloorani et al., 2008;Julien & Sobrino, 2019;Kandasamy et al., 2013;Moreno et al., 2014;Sarafanov et al., 2020;Siabi et al., 2020;Weiss et al., 2014) or remotely sensed soil moisture (Wang et al., 2012).However, the above-mentioned gap-filling procedures have disadvantages.The implicit gap-filling treatment creates regional inconsistencies, as the resulting monthly maps are obtained with inconsistent observation quantities and timings.Furthermore, it distorts the missing value pattern of the data set.Both explicit and implicit missing value treatments are usually univariate, which means that they rely exclusively on information from the target variable.These approaches thus cannot account for information from other variables that might have been observed simultaneously in the same area.Finally, the complex missing value pattern of satellite-retrieved variables depends on measurement limitations, which depend on the observed data, a feature called "missing not at random" or MNAR (Little & Rubin, 2014;Rubin, 1976).For example, satellite-observed LST derived via thermal infrared emission of the land surface cannot measure temperature below clouds.However, LSTs below clouds can be systematically warmer or colder than the observed temperatures under clear-sky conditions in their immediate neighborhood.Accounting for this unknown bias in incomplete data with MNAR missingness is difficult with simple gap-filling procedures.
The literature on multivariate gap-filling methods applied to Earth observations is sparse.Moreno-Martínez et al. (2020) gap-fill and fuse remotely sensed MODIS and LANDSAT optical bands by applying a pixel-wise linear regression model between both data sets.Kadow et al. (2020) gap-fill historical temperature observations by searching similar conditions in the modeled record ("analog search") via a neural network approach attempting image inpainting.K. Liu et al. (2023) gap-fill ESA-CCI surface layer soil moisture with a sophisticated neighborhood search algorithm and information from land surface models.Several studies gap-fill remotely sensed soil moisture by employing machine learning approaches that learn the relationship between observed soil moisture and related variables, for example, temperature, precipitation, Normalized Difference Vegetation Index, or

Data
Within this study, we demonstrate a use case of CLIMFILL using real-world observations.The goal is to collect a set of satellite and ground observations, gap-fill them using CLIMFILL, and publish the resulting data set.For that we consider two different sets of data: first, observational products are included as an input for the gap-filling framework that will comprise the end product.Second, independent data sets are used to benchmark the gap-filling.

Input Data
The input data consist of global, gridded, monthly data products representing relevant land surface climate characteristics.These can be from satellite remote sensing or ground observations.The included data sets are listed in Table 1.We focus on the time period 1995-2020.Data sets are regridded to a common 0.5° degree grid.
We furthermore include gridded observations of precipitation and temperature to simulate ground observations.Ground observations are usually independent of satellite observations and offer 2-m temperature, which is not observable from space.Sparsely located ground observations are inherently similar to gridded observations with many missing values.For simulating ground observations, we use CRU TS precipitation and 2-m temperature (I.Harris et al., 2020).We take the gridded data set and remove all grid points that do not contain a station.In doing so, we circumvent the problem of representativeness between point-scale ground observations and gridded satellite observations and can still assess the merit of the multivariate gap-filling to station data.Since CRU station locations are only available until version 3.26, which extends until 2017, we elongated the station locations from 2017 until the end of 2020, assuming that no stations were abandoned or added between 2017 and 2020.

Benchmarking Data Sets
Naturally, the gap-filled data set cannot be compared to the original observations, as the original values are missing.To evaluate the overall characteristics and benchmark the gap-filled data set, we, therefore, compare the estimates of CLIMFILL to the corresponding variables of the ERA5-Land global reanalysis.ERA5-Land is a state-of-the-art reanalysis product providing gap-free maps of many relevant surface variables.It is produced by running the land model of the European Centre for Medium-Range Weather Forecasts (ECMWF) with downscaled atmospheric forcing from the ERA-5 reanalysis, to create a land-only reanalysis with higher resolution and in better agreement with observations of land hydrology, such as soil moisture (Muñoz-Sabater et al., 2021).Of all the variables considered for gap-filling (Table 1), terrestrial water storage (TWS; and groundwater) and BA are not available in ERA5-Land.Some modeling approaches exist that tweak atmospheric or land models to produce estimates of groundwater storage (Döll et al., 2014;Hanasaki et al., 2018;Pokhrel et al., 2015), but this is not yet implemented in ERA5-Land (Muñoz-Sabater et al., 2021).ERA5-Land does not have a fire module; therefore, BA is also not available and hence cannot be benchmarked.None of the input data described in Table 1 is directly dependent on ERA5-Land or similar reanalysis products.Note also that precipitation estimates from ERA5 have been shown to diverge substantially from other observation-based sources (Alexander et al., 2020) and therefore should be used only carefully.Furthermore, gap-filled ESA-CCI surface layer soil moisture is compared to ground observations of soil moisture included in the International Soil Moisture Network (ISMN; W. Dorigo et al., 2021;W. A. Dorigo et al., 2011).The ISMN is the largest volunteer-based collection of measurements of soil moisture stations across the globe that has been homogenized, quality checked, and freely available online (W. A. Dorigo et al., 2011).

Infilling Missing Values Using CLIMFILL
CLIMFILL is a multivariate gap-filling framework that combines state-of-the-art spatial interpolation with an iterative approach accounting for the dependencies across multiple incomplete variables.CLIMFILL is described and benchmarked in detail in Bessenbacher et al. (2022), but we give a short overview here and also explain updates to the framework employed in this study.
The CLIMFILL workflow is divided into four steps.In the first step (initial guess), all gaps are filled by spatial interpolation, creating an initial guess of the missing values.This interpolation follows the procedure from Haylock et al. (2008), where the monthly climatology is gap-filled using thin-plate-spline interpolation, and the anomalies are gap-filled using a kriging procedure that is optimized for large data.In the second step (feature engineering), lagged features from the available variables are created additionally.They comprise 3-month, 6-month, and 12-month running averages both backward and forward in time, to include information about slow-changing

1901-2017
Grid points without stations I. Harris et al. (2020) Note.Variables with an asterisk (*) have been subject to monthly aggregation, where monthly means are computed values if at least 50% of the days within this month contain a valid observation (for more details see text).Variables marked with a circle (°) are ground observations, and all others are satellite observations.

Table 1
Remotely Sensed, Incomplete Satellite Observations and Sparse Ground Observations Used for Multivariate Gap-Filling processes like seasonality or soil moisture storage.In the third step (clustering), the data points are divided into clusters using the K-Means algorithm.This happens independently of the temporal and spatial position of the values and only accounts for their multivariate similarity.In the last step, the multivariate context of each missing value is used to update the initial guess from interpolation.In this fourth step (regression step), the initial estimates of step 1 for the missing values are iteratively updated in each cluster using an adapted version of the MissForest algorithm (Stekhoven & Bühlmann, 2012) with a Random Forest Regressor until the estimates for the missing values do not change anymore between iterations (convergence) or a maximum number of iterations is reached (early stopping).
For this study, we have adjusted the initial gap-fill to better account for the individual missingness patterns of the different variables.The following procedure is applied to obtain the initial guess of step 1: 1.For variables with a daily resolution, but with many missing values (see Table 1; marked with an asterisk), we use the monthly mean of all available measured values as the initial guess, if less than 50% of the days are observed within this month.For ESA-CCI BA, only the fraction of the observed area and not the fraction of observed days per month is provided.We apply an equivalent criterion, where observations are valid if at least 50% of the area is observed.Diurnal temperature range (DTR) is computed as the difference between MODIS day and night data, wherever both are available.If more than 50% of the days are observed, the value is not set to missing and therefore not gap-filled.2. For variables with a shorter time coverage than 1995-2020 or whole months missing, the initial guess for the missing months and years is the mean seasonal cycle of all observed years.3.All remaining missing values are filled using the spatial interpolation procedure as described in detail in Bessenbacher et al. (2022).

Validation and Cross-Validation
Since we do not know the original value of the gaps in the data sets, we need a procedure to validate our estimates and for optimizing the hyperparameters of the interpolation and regression steps described above.Therefore, originally observed data are deleted in the years 2004 and 2005 for each variable.These values are gap-filled with CLIMFILL and used for hyperparameter estimation and for assessing the performance of the gap-filling.
Since randomly missing points are relatively easy to gap-fill through spatial or temporal interpolation, we delete data according to a pattern of randomly placed rectangles of 18° latitude and 9° longitude (later referred to as "validation rectangles," Figure 1) that are differently placed for each variable and month.We create 10 sets of validation rectangles, such that in each of the validation sets approximately 10% of the data are missing within the validation rectangles.These 10 sets of validation data are gap-filled individually and combined in the end such that all observed data points within the time period are missing and gap-filled in one of the 10 sets.With this 10.1029/2023JD039099 6 of 19 procedure, we ensure that the created gaps emulate real gaps in satellite data in two properties, namely that multivariate dependencies can still be learned because the missingness patterns are not congruent and spatial interpolation is difficult because the gaps are large enough such that some grid points are far away from the spatially closest observed point.All validation sets are used for the validation of CLIMFILL in Section 4.2, and one of the sets is used for cross-validating the hyperparameters to avoid an unmanageable computational expense.The cross-validation of the initial guess and the regression step, that is, for hyperparameter settings of the thin-platespline interpolation, the kriging, and the random forest is done by repeatedly gap-filling the validation rectangles in one of the validation sets using a different combination of hyperparameters each time.The gap-filled data are then compared to the originally observed one and the hyperparameter combination with the lowest root mean square error (RMSE) in the gap-filling is chosen.

Missing Value Patterns
The spatiotemporal distribution of missing values in the included data sets is shown in Figure 2. ESA-CCI surface layer soil moisture (Figure 2a) has a high amount of missing values: it is missing in approximately 59% of the cases.Missing values occur especially in the tropics, where dense vegetation leads to the signal being dismissed in postprocessing since microwave remote sensing rather measures canopy water content than soil water content.Furthermore, retrieval is not possible when there is snow cover, such as in high latitudes in winter.Surface temperature from MODIS (Figure 2b) cannot be observed below clouds, which leads to high fractions of missing values in the inner tropics and the extratropical storm tracks.Surface temperature is only observed in approximately 42% of the cases, making it a similarly challenging case for gap-filling as soil moisture.Precipitation from GPM (Figure 2c) is missing only in high latitudes in winter.TWS from GRACE (Figure 2d) is missing if individual months of measurements do not converge to a valid signal in the retrieval algorithm.Furthermore, there is a gap in 2017 between the two GRACE and GRACE-FO missions.ESA-CCI snow cover fraction (SCF; Figure 2g) cannot be observed below clouds and in the polar night, leading to 68% of missing values.DTR (Figure 2h) is missing if both night and day measurement of the MODIS instrument is hampered by clouds, leading to the highest amount of missing values among all remotely sensed observations (68%).Analog to surface temperature, BA (Figure 2i) is only retrievable under clear-sky conditions, which leads to approximately 59% of the values missing, especially in areas that have a significant fire regime.
As mentioned above, CRU temperature and precipitation (Figures 2e and 2f) are observed in grid points that contain an active CRU station.Those generally decrease over time and there are fewer stations included measuring temperature than precipitation.Counting grid points that contain an active station in time and space, ground observations of temperature and precipitation have the highest overall fractions of missing values (92% and 84%, respectively).
For multivariate analysis, satellite data that are not gap-filled can only be used if all considered observations are available at the region and time of interest.In Table 2, we show the fraction of data points for any combination of two variables where both are observed at the same time.For an exemplary analysis of surface temperature and soil moisture dynamics, only 21% of the data points are available, which makes the analysis difficult.Due to the missingness patterns of both variables, only clear-sky vegetation sparse snow-free cases can be investigated, which do not cover all dynamics of land-atmosphere interaction of interest.This is similar for other combinations of variables, as only TWS in combination with remotely sensed precipitation or BA features a larger cooccurrence fraction than 50%.These values naturally decline even further if more than two variables would be taken into account.Figure 3 shows the number of observations available at each grid point each month, averaged over the 1996-2020 timespan.It shows that in many regions of the world, only a small subset of variables are available at the same time, and the inner tropics are especially underobserved.This renders multivariate analysis difficult.

CLIMFILL Validation
The CLIMFILL gap-fill is validated by deleting and subsequently gap-filling originally observed data in the "validation rectangles" (see Section 3.2).The performance of the gap-fill is compared with the first step of CLIMFILL, the initial guess, to show the added value of the multivariate gap-filling as compared to univariate interpolation.Figure 4 shows the Pearson correlation and RMSE of the initial guess and the CLIMFILL estimates with the original observations.Additionally, Figures 5 and 6 show the RMSE and the Pearson correlation coefficient on anomalies for the CLIMFILL gap-fill averaged over the IPCC regions.As there are not always observations available for each grid point, some regions cannot be validated.In the following, we will discuss the merit of the CLIMFILL gap-fill by referring to the overall performance shown in Figure 4 and its spatial structure of error and correlation in Figures 5 and 6.
For both correlations and RMSE, the gap-fill estimates for ESA-CCI surface layer soil moisture (SM) improve for the whole time series, the anomalies, and the seasonal cycle, compared to the initial guess interpolation (Figure 4).This improvement is most pronounced for the anomaly correlation and has no pronounced spatial pattern (Figures 5a and 6a).Correlation on absolute values is higher compared to the spatiotemporal ESA-CCI surface layer soil moisture in Llamas et al. (2020) and comparable to the spatiotemporal gap-fill of Wang et al. (2012), although they delete validation points randomly, which is significantly easier to gap-fill than larger spatiotemporal gaps like the validation rectangles used in this study (Bessenbacher et al., 2022).MODIS LST (Figures 4, 5, and 6b) and MODIS DTR (Figures 4, 5, and 6h) gap-fill improve compared to the first guess only for the anomalies.This generally indicates that the spatial structure and the strong seasonality of these variables favor interpolation.However, the added value of the multivariate gap-filling with CLIMFILL for the monthly anomalies shows that at the scale of weather dynamics, multivariate dependencies can inform gap-filling of ground temperature well.This might also be visible in the higher correlations in the tropics, where seasonality is less pronounced.The seasonal change of MODIS temperature is quite smooth and therefore spatiotemporal interpolation is already quite a good guess.
Gap-fill estimates of GPM precipitation (PSAT; Figure 4) also improve for all scores on all-time series, in particular for the anomaly time series with no distinct spatial pattern (Figures 5 and 6c).Note that a previous study, this variable was the most challenging case for gap-filling (Bessenbacher et al., 2022).Bessenbacher et al. (2022) argued this is due to the high local variability and logarithmic distribution of precipitation-here, this might be alleviated by the fact that monthly precipitation is considered and only originally observed values are deleted, such that gap-filling GPM precipitation (for which high-latitude winter gaps are not filled) becomes easier.PSAT gap-fill generally performs worse in high latitudes, where the missingness is higher (compare Figure 2) and in tropical South America.

GRACE TWS (Figures 4, 5, and 6d
) is only missing in specific months, where the seasonal cycle gap-fill in the first step of CLIMFILL is already quite a good guess for overall patterns.Interestingly, its gap-fill has a distinct regional pattern and the highest error in the tropics, although the fraction of missing values is uniform across the globe.4, 5, and 6f) have a lower correlation and higher RMSE than the initial guess in most cases.Those station observations prove to be a difficult case for multivariate gap-filling: both because their missingness is high, and spatial interpolation is already quite a good guess.Their regional errors correspond to the station distributions: as there are little available observations in the tropics, the gap-fill errors are higher and the correlations are lower.
ESA-CCI SCF (Figures 4, 5, and 6g) has better performance than the initial guess for RMSE but not for correlations.Both correlation and RMSE have high regional variability.This could be due to the fact that ESA-CCI SCF contains a lot of zeros, making interpolation easy and creating spurious relations within the multivariate gap-filling step.
Overall, ESA-CCI BA has both the weakest correlations and the highest RMSE for both the initial guess and the CLIMFILL estimate (Figure 4), highlighting that this variable is influenced by not only climate processes but also ecosystem dynamics and human management.Its errors are highest in regions with a strong fire regime.In particular, BA depends on the other variables via a complex multiscale process including fuel buildup, vegetation biomass, ignition, and fire management; it does not show a clear deterministic response to the atmospheric forcing described by the other variables, which might impact the gap-filling.Interestingly, the CLIMFILL estimate has not only a lower RMSE but also lower correlation than the initial guess.In summary, the effectiveness of CLIMFILL's gap-filling approach varies significantly between the variables and regions under consideration.In particular, CLIMFILL outperforms univariate interpolation especially for ESA-CCI surface layer soil moisture (SM), GPM precipitation (PSAT), and GRACE TWS for the majority of scores.This shows that the multivariate gap-filling step is valuable in these cases.The largest merit is shown in anomalies, where weather dynamics play a substantial role and multivariate relationships are strong.The situation differs for MODIS LST, MODIS DTR, ESA-CCI BA, and ESA-CCI SCF, where the results are mixed, highly dependent on the score, the region, and/or the time aggregation in question.This case can still provide useful results but needs to be assessed with caution.Finally, for station-based CRU temperature and precipitation data, the multivariate gap-filling approach does not yield significant advantages over the univariate initial guess.Nevertheless, these variables can still play a supportive role in enhancing the gap-fill process for other included variables.
We conclude that the effectiveness of CLIMFILL's gap-filling method is context-specific, with its success dependent on the nature of the gaps and the variable.However it has shown to have superior results compared to interpolation for variables with high fractions of missing values and complicated missingness patterns like surface layer soil moisture, precipitation, and TWS.

Comparison to ERA5-Land Reanalysis
We compare the original, incomplete data sets and the data sets gap-filled with CLIMFILL to ERA5-Land reanalysis data.This serves to show the impact of ignoring gaps in observations on the representation of anomalies and trends within the observed data.To demonstrate this impact, we create averages of the variables aggregated over regions.These regional averages consist only of the observed values in the original, incomplete data set, but are gap-free in the CLIMFILL product and ERA5-Land.Therefore, the averages consist of a different number of observations for the incomplete original data.In this way, the impact of ignoring gaps while creating spatially aggregated estimates of anomalies and trends can be shown.Figures 7a-7i show maps of the correlation between the time series of normalized anomalies between ERA5-Land and the CLIMFILL estimates calculated for each IPCC region.Generally, correlations are high (>0.8)for most variables.Lower correlations for the variables generally occur where the fraction of missing values is high (compare Figure 2).ESA-CCI SCF correlations are small in tropical regions, but in those regions, it can mostly be assumed to be zero.Unfortunately, no estimates of TWS or BA exist in ERA5-Land.
Figures 7j and 7k show the Pearson correlation coefficient and RMSE for both the CLIMFILL estimates and the original, incomplete data sets compared to ERA5-Land.As mentioned above, the scores are calculated by first averaging all data sets over the IPCC regions on land, then applying the scores on the aggregated time series, and finally creating boxplots of the scores across IPCC regions.CLIMFILL estimates have a higher correlation and lower RMSE with ERA5-Land for all variables except GPM precipitation and CRU 2-m temperature, showing gap-filling with CLIMFILL in most cases is beneficial for representing anomalies in observations if compared to ERA5-Land.ERA-Interim precipitation has shown to diverge from observations quite substantially (Alexander et al., 2020), but Hersbach et al. (2020) report increased agreement from ERA-Interim to ERA5 to the Tropical Rainfall Measuring Mission (TRMM), which is included in the GPM data set (Huffmann et al., 2019).Therefore we conclude that the gap-fill for precipitation has a lower agreement with ERA5-Land as its high spatial variability and the skewed distribution make gap-filling quite challenging.CRU 2-m temperature gap-filling is solely happening outside of existing stations and might be hampered by the low station density, which leads to the highest fraction of missing values among all considered data sets.

Comparison to ISMN Soil Moisture Data
We further compare original and gap-filled ESA-CCI surface layer soil moisture to independent ground observations of soil moisture.The ISMN collects observations from in situ soil moisture stations around the world (W.Dorigo et al., 2021;W. A. Dorigo et al., 2011).Figures 9a and 9b show the correlation of original and gap-filled ESA-CCI surface layer soil moisture observations with the soil moisture measured at the ISMN stations at the respective grid point.The correlations are similar, despite a bit lower for the gap-filling compared to the original observations.Other studies attempting gap-filling of ESA-CCI surface layer soil moisture find similar agreement of incomplete and gap-filled remotely sensed soil moisture with ISMN station observations (Almendra-Martín et al., 2021;Y. Liu, Yang, et al., 2020;Llamas et al., 2020).Figures 9c-9e show the time series of soil moisture observations at the two ISMN stations with the longest data record, and at the ISMN station with the longest record outside the USA.
Despite the overall slightly worse performance in individual locations (cf., Figures 9a and 9b), overall, the correlations of data gap-filled with CLIMFILL compared to interpolated data (initial guess) are higher in all months and this difference is most pronounced in winter, when data are scarcest (Figure 10).Furthermore, gaps in ESA-CCI might occur in more difficult measuring conditions (e.g., winter, dense vegetation) which could impact the correlations of the gap-filled data versus the in situ data.But also the quality of the in situ may be impacted by for example, variations in soil temperature (see e.g., Mittelbach et al., 2012).We finally argue that the main benefit of CLIMFILL is its ability to plausibly elongate the ESA-CCI time series (and other observations).By filling the observational record this way, the number of datapoints where these different observations can be observed is increased massively 10.1029/2023JD039099 13 of 19 from around 10% to 30% depending on the season to 100% of all station observations (Figure 10).This allows for more robust analyses such as satellite-based trend estimates (e.g., Golian et al., 2019;Zhao et al., 2019).

Illustrative Example: Major Fire Seasons Around the Globe
CLIMFILL stands out as a unique tool capable of delivering gap-filled data within the realm of multivariate analysis.The capabilities of CLIMFILL extend beyond the boundaries of this current study, offering the potential for a plethora of future investigations.One exemplary illustration of its use arises within the realm of compound event analysis, a domain where diverse data sources of many variables must be seamlessly harmonized, for example, through multivariate gap-filling techniques.In Figure 11, we show the analysis of the record-breaking, catastrophic fire season 2017 in Europe (European Commission Joint Research Centre, 2022).We show gap-filled BA (Figure 11c) and the percentiles of heat and moisture-related variables (Figures 11d-11j), effectively using almost all the CLIMFILL estimates variables.Temperature anomalies in the uppermost percentiles (0.9 or higher) or unprecedented temperatures are visible for all events at 2 m and at the surface.The DTR is also above-average high.Severely desiccated soils from the lowest percentile precipitation amounts are also visible in the data.A major advantage of the CLIMFILL estimates is that the gap-filling is multivariate, and this leads to a better representation of anomalies outside the normal seasonal cycle and a higher, spatially smooth field of correlations between soil moisture and temperature, showing increased representation of land-atmosphere interactions (Figures 11a and 11b).

Discussion and Conclusions
Combining different remotely sensed and in situ observations to create a coherent view of the land climate is an active field of research and essential for the multivariate multisource analysis of land-atmosphere coupling and compound events.Large fractions of missing values and complex patterns of missing values impede such endeavors as the numbers, sizes, and locations of the gaps are often not overlapping.Implicit gap-filling attempts conflate different numbers and times of observations into one monthly grid, and explicit gap-filling treatments usually cannot take into account cross-variable dependency or are not independent of physical models.We apply the recently developed spatiotemporal and multivariate gap-filling algorithm CLIMFILL to nine remotely sensed and in situ observations centering on land-climate processes, especially heatwaves and droughts.CLIMFILL leverages temporal autocorrelation, spatial neighborhood, and cross-variable dependence to estimate missing values in multivariate Earth system data.By doing so, we created a spatially complete data set, consisting of a range of observations, without strong physical model assumptions inherent to reanalysis data.
The considered observations have a large fraction of missing values (22%-92%, depending on the variable) and low overlap in their missingness patterns, leading to an overall minuscule fraction of grid points observed by several variables at the same time.Additionally deleted, observed points serve as an out-of-sample validation for the gap-filled data set, showing an overall clear edge of the CLIMFILL algorithm compared to state-of-theart interpolation, especially when looking at monthly anomalies.This highlights the added value of considering cross-variable dependence and thus a larger fraction of the available information when estimating dynamical properties of unobserved Earth system dynamics.
We benchmark the merit of the gap-filling with other, independent data sets.Comparing the gap-filled observations with ERA5-Reanalysis shows high correlations and low RMSE on monthly anomalies and higher agreement with ERA5-Land as the original, incomplete data sets for most variables.Furthermore, trends in ESA-CCI surface layer soil moisture only agree with ERA5-Land trends for the Northern Extratropics after gap-filling, and the wetting trend of the incomplete data set can be attributed to a decreasing number of missing values over time.These results indicate that ignoring gaps or implicit gap-filling introduces biases and can change trend signs compared to ERA5-Land which are alleviated with multivariate gap-filling.Comparison of ESA-CCI surface layer soil moisture with station observations from the ISMN shows overall high correlations and can help elongate remotely sensed observations to ease the comparison with in situ stations and facilitate other analyses.A multivariate analysis of major global fire events shows a consistent picture with strong anomalies in dryness and temperature for the considered variables.The CLIMFILL algorithm thus allows recovery of a gap-free picture of these extreme events, whose analysis is usually hampered by incomplete data.CLIMFILL is able to mature further in the future by including more variables and increasing the temporal and spatial resolution.
Overcoming the fragmented, incomplete record at the land-atmosphere interface is a crucial step for multivariate, multisource analysis, for example, of extreme or compound events like heatwaves, droughts, and wildfires.Ignoring gaps in observations when aggregating data can introduce biases to the observational record and prevent multivariate analysis.Multivariate data-driven gap-filling at the heart of such analysis can improve the representation of trends, anomalies, and multivariate dependency structure of independent, originally incomplete observations.In the buildup toward a Digital Twin of the Earth, multivariate gap-filling can take a central role.

Figure 1 .
Figure 1.Illustration of the validation procedure.Monthly mean surface layer soil moisture for September 2005 from European Space Agency (ESA)-Climate Change Initiative (CCI) soil moisture with gaps from measurement limitations (unobserved points) and additionally added gaps for validation based on "validation rectangles" (validation points).

Figure 2 .
Figure 2. Fraction of missing values per variable in the observations.Maps show the fraction of missing values per grid point on land, rectangular plots show the time evolution of the fraction of missing values per latitude (y-axis) and month (x-axis).

Figure 4 .
Figure 4. Validation scores of the gap-filling.Boxplots showing the spatial distribution of Pearson correlation and normalized root mean square error (RMSE) between validation points of the initial guess and CLIMate data gapFILL (CLIMFILL) estimates for all variables.The scores are calculated on the original time series, the anomalies, and the mean seasonal cycle separately for each grid point and then summarized into boxplots.All variables are scaled to mean zero and standard deviation one.The boxplots show the median as a black horizontal line, the boxes are the quartiles, and the whiskers show 1.5 times the quartile length.

Figure 5 .
Figure 5. Pearson correlation of the gap-filling.Maps of Pearson correlation between the originally available observations within the validation rectangles and the CLIMate data gapFILL (CLIMFILL) gap-fill.Correlations are calculated with the anomalies, that is, the mean seasonal cycle is removed.

Figure 8
Figure 8 shows linear trends in soil moisture in (a) the original ESA-CCI surface layer soil moisture, (b) the gap-filled ESA-CCI surface layer soil moisture and in ERA5-Land, (d) both gap-free and (c) the part that is observable in the ESA-CCI data set (denoted "gaps deleted").Trends in ERA5-Land are generally larger than in ESA-CCI.Gap-filling ESA-CCI surface layer soil moisture leads to less noisy trends, especially in regions

Figure 6 .
Figure 6.Root mean square error (RMSE) of the gap-filling.Maps of normalized RMSE between the originally available observations within the validation rectangles and the CLIMate data gapFILL (CLIMFILL) gap-fill.Correlations are calculated with the anomalies, that is, the mean seasonal cycle is removed.The values are normalized (mean zero, standard deviation one) before RMSE calculation.

Figure 7 .
Figure 7. Comparing the CLIMate data gapFILL (CLIMFILL) gap-filled data and the incomplete original data to the ERA5-Land reanalysis: (a-i) Pearson correlation of anomalies of all gap-filled data sets between ERA5-Land and the CLIMFILL estimates aggregated over IPCC regions.(j) Pearson correlation and (k) root mean square error (RMSE) compared to ERA5-Land.The scores are calculated on the normalized anomalies, that is, when the seasonal cycle is removed and all variables are scaled to mean zero and standard deviation one.Scores are calculated separately for each IPCC region, the bars correspond to the median of all IPCC regions, and the uncertainty ranges are the standard deviation across the regions.

Figure 8 .
Figure 8. Trends in soil moisture for (a) original and (b) gap-filled European Space Agency (ESA)-Climate Change Initiative (CCI) surface layer soil moisture as well as (d) original ERA5-Land surface layer soil moisture and the same without the data points that are not originally observed in ESA-CCI (c).Trends in maps (a)-(d) are calculated for each grid point as the slope of a linear trend in the years 2009-2020 on yearly averages.(e) Yearly time series averaged over the Northern Extratropics.

Figure 9 .
Figure 9.Comparison with in situ soil moisture observations.Correlation between International Soil Moisture Network (ISMN) in situ soil moisture observations and European Space Agency (ESA)-Climate Change Initiative (CCI) surface layer soil moisture at the grid points containing the station globally.The correlations are computed on the subsets of time series at the stations that are observed (a, correlation between gray and green lines in (c)-(e)) or gap-filled (b, correlation between orange and green lines in (c)-(e)), that is, for each station the correlations are computed on two nonoverlapping time series.(c-e) Time series of soil moisture at three sample stations with a long observational record in the ISMN, together with incomplete and gap-filled ESA-CCI surface layer soil moisture at the respective grid point.Displayed are all ISMN stations that have at least 5 years of data.

Figure 10 .
Figure 10.Impact of missing values on correlations.Average number of observed months in the International Soil Moisture Network (ISMN) data from 1995 to 2020 (gray bars) and mean correlation between the ISMN data and the European Space Agency (ESA)-Climate Change Initiative (CCI) surface layer soil moisture gap-filled by interpolation (blue) and CLIMate data gapFILL (CLIMFILL; orange).Displayed are all ISMN stations that have at least 5 years of data.

Figure 11 .
Figure 11.Multivariate analysis of a major European fire season in August 2017.Representation of event in original ("With Gaps") and gap-filled data ("Gap-Filled") for (c) burned area and (d-k) percentiles of the respective variable in dryness/hotness reached in August 2017.A value of one means the conditions are unprecedented at this location.(a, b) The correlation between surface layer soil moisture and surface temperature for all August months between 1995 and 2020, for original and gap-filled data, respectively.

Table 2
Cooccurrence of All Combinations of Two Observations Globally in the Time Period 1995-2020, Expressed as the Percentage of Data Points Where Both Observations Are Available Figure 3. Cooccurrence of variables.The average number of available monthly observations per grid point for the 1995-2020 time period.For example, areas with the value 4 have, on average, four observations per month that are available, while the rest (5) are missing.