Synergetic retrieval from multi-mission spaceborne measurements for enhanced aerosol and surface characterization

Atmospheric aerosol is one of the main drivers of climate change. Currently, a number of different satellites in Earth orbit are dedicated to aerosol studies. Due to limited information content, the primary aerosol product of most satellite missions is AOD (Aerosol Optical Depth), while the accuracy of aerosol size and type retrieval from spaceborne remote sensing still requires improvement. Combining measurements from different satellites increases their information content and, therefore, can provide new possibilities for retrieving an extended set of both aerosol and surface properties.

In this paper, we present the physical basis and concept of the recently developed synergetic approach for aerosol and surface characterization using diverse spaceborne measurements (hereinafter SYREMIS (SYnergetic REtrieval from Multi-MISsion instruments) approach). The approach was implemented in the GRASP (Generalized Retrieval of Atmosphere and Surface Properties) algorithm and has been tested on two types of synergetic measurements: (i) synergy of polar-orbiting satellites (LEO + LEO synergy combining Sentinel-5P/TROPOMI, Sentinel-3A/OLCI, and Sentinel-3B/OLCI instruments), (ii) synergy of polar-orbiting and geostationary satellites (LEO + GEO synergy based on Sentinel-5P/TROPOMI, Sentinel-3A/OLCI, Sentinel-3B/OLCI, and Himawari-8/AHI instruments). On the one hand, such a synergetic satellite constellation extends the spectral range of the measurements. On the other hand, it provides unprecedented global spatial coverage with high temporal resolution, which is crucial for a number of climate studies. It is shown that the SYREMIS/GRASP approach facilitates the transfer of information content from instruments with richer information content to those with lower one. This results in substantial enhancements in aerosol and surface characterization for all instruments within the synergy.

1 Introduction

Many climate studies require global extended aerosol and surface characteristics, including properties such as Aerosol Optical Depth (AOD) and size distribution, Single Scattering Albedo (SSA), full surface Bidirectional Reflectance Distribution Function (BRDF), etc. This is particularly relevant for the generation of high-quality aerosol and surface Essential Climate Variables (ECVs), which characterize the Earth's climate system, as well as for air-quality monitoring, aerosol emission and transport studies (Dubovik et al., 2008, 2021b; Martin, 2008; Hollmann et al., 2013; IPCC, 2021; Chen et al., 2019, 2022a). In addition to a global scale, the high temporal resolution of the extended aerosol properties is required for such important but challenging studies as aerosol-cloud interactions, gas-to-particle transformation, and atmospheric aerosol dynamics (Pöschl, 2005; Rosenfeld et al., 2023; Vehkamäki and Riipinen, 2012).

Global information about aerosols can be obtained from spaceborne measurements. Therefore, climate studies increasingly rely on high-quality aerosol characterization from space. At present, there are many different satellites in Earth orbit dedicated to aerosol characterization. Nevertheless, due to their limited information content, the primary aerosol product of most satellite missions is AOD (Remer et al., 2005; Levy et al., 2013; Sayer et al., 2018b; Sogacheva et al., 2020), while the accuracy of extended properties with high temporal resolution still requires improvement. To address this problem, several requirements for satellite measurements have been formulated based on the general principles of light scattering theory and atmospheric dynamics studies (Van der Hulst, 1981; Tsang et al., 1985; Bohren and Huffman, 1998; King et al., 1999; Mishchenko et al., 2002, 2004; Hasekamp and Landgraf, 2007; Hasekamp et al., 2011, 2024; Lenoble et al., 2013; Dubovik et al., 2011, 2019; Remer et al., 2024). In particular, such measurements should include:

• i.

  Multi-angular measurements in a wide range of scattering angles where the differences between the angular dependence of aerosol and surface signals can be observed, and angular sampling is sufficient for aerosol characterization.

• ii.

  Measurements in a wide spectral range (preferably from the Ultraviolet (UV) to Shortwave Infrared (SWIR) ranges) to take advantage of the different spectral dependence of aerosol and surface signals and to observe spectral features of different aerosol species.

• iii.

  Polarimetric measurements to exploit differences in the polarization signatures of aerosol and surface signals and the strong dependence of such measurements on microphysical properties of aerosols.

• iv.

  Frequent temporal measurements to account for the temporal variability of aerosol properties, as well as differences in aerosol and surface conditions.

Conditions (i)–(iii) are substantially covered by Multi-Angular Polarimetric (MAP) spaceborne missions such as, for example, POLDER-3/PARASOL mission, which ended in 2013 (Deschamps et al., 1994; Tanré et al., 2011), or new missions, like PACE (with HARP-2 and SPEX instruments on board), MetOp-SG 3MI, and CO2M MAP (Dubovik et al., 2019; Hasekamp et al., 2019; McBride et al., 2024; Werdell et al., 2024; Fu et al., 2025; Sienkiewicz et al., 2025). Nevertheless, each of these spaceborne instruments, including MAP, still has limitations related to their spectral/spatial coverage and/or resolution. Moreover, most polar-orbiting satellites have a relatively low revisiting frequency (a few times per day or less), which is insufficient for studying physical, chemical, and dynamic processes in the atmosphere with the required temporal resolution of one hour or even better.

Strictly speaking, none of the currently operating and future aerosol-oriented satellite instruments alone completely meet all the aforementioned requirements (i)–(iv). The solution to this problem has been discussed for a long time and is based on the synergetic retrieval of combined measurements from different sensors (Aires et al., 2012; Holzer-Popp et al., 2008; Vanhellemont et al., 2014; Wang et al., 2014; Lee and Ahn, 2021). Indeed, different satellites dedicated to atmospheric studies may have varying spectral coverage and resolution, observe the same area on Earth's surface on the same day but at different times or relative positions. As a result, when properly collocated and combined, such measurements can provide multi-angular, multi-temporal measurements in an extended spectral range, satisfying requirements (i)–(iv).

Despite the known approach, a generalized synergetic retrieval method applicable to diverse multi-instrument L1 measurements (L1 synergy) has not yet been developed. The problem is algorithmic rather than engineering. To exploit the combined L1 synergetic measurements, the retrieval algorithm must satisfy the following conditions:

• v.

  The retrieval should be based on flexible forward models adaptable to the information content of the measurements.

• vi.

  The retrieval should be able to account for diverse measurements with potentially different calibration accuracy, spectral, and spatial resolutions.

• vii.

  The algorithm should be able to account for multi-temporal (not collocated in time) measurements.

Many retrieval algorithms can meet requirements (v) and (vi), resulting in synergetic approaches for collocated-in-time satellite measurements. Examples include the synergy of MERIS and AATSR from the ENVISAT platform (North et al., 2008), the synergy of OLCI and SLSTR from the Sentinel-3 platform (Henocq et al., 2018), and the PMAP synergetic algorithm for GOME-2, AVHRR, and IASI on the MetOp platform (Grzegorski et al., 2021). Nevertheless, the correct treatment of observations that are not collocated in time is still beyond the capacity of most existing algorithms. Since aerosol properties do not change randomly in time and space, showing temporal and spatial correlations due to atmospheric dynamic processes, accounting for such temporal dependencies is crucial in synergetic retrieval.

Here, we use the GRASP (Generalized Retrieval of Atmosphere and Surface Properties) algorithm (Dubovik et al., 2011, 2014, 2021a), which has emerged from the successful heritage of AERONET retrieval developments (Dubovik and King, 2000; Dubovik et al., 2000, 2002; Dubovik, 2004) and pursues the idea of creating a scientifically rigorous and versatile algorithm. As a result, GRASP has significantly extended capabilities and areas of applicability. In particular, it can be applied to diverse remote sensing observations (passive, active, ground-based, satellite, in situ, etc.), can simultaneously retrieve a large number of different atmospheric characteristics, and is well-optimized for synergetic retrievals. This has been achieved by pursuing the following generalization principles:

• The two main modules of the algorithm, the “Forward model” and “Numerical Inversion,” are independent. This allows for continuous development and extension of all functionalities of both modules without compromising previously established applications.

• The “Forward model” provides consistent modeling of all observations to which GRASP can potentially be applied.

• The “Numerical Inversion” is general and flexible enough for inverting different observations and retrieving all atmospheric and surface parameters that affect those observations.

GRASP relies on a statistical optimization approach based on the Multi-Term LSM (Least-Square-Method) (e.g., Dubovik et al., 2021a). This approach, in contrast to the more common Optimal Estimation (Rodgers, 2000), emphasizes the use of multiple a priori constraints. Indeed, the retrieval of most atmospheric parameters from remote sensing observations is an ill-posed problem, and the use of a priori constraints is necessary for successful retrieval. However, such constraints can be very different for various retrieved atmospheric characteristics (e.g., aerosol size distribution, vertical profile, index of refraction, surface reflectance parameters). A single constraining approach (such as the direct use of a priori estimates for each retrieved parameter) is not optimal and is hardly possible in some situations. Using Multi-Term LSM allows resolution of this difficulty by using different a priori constraints for diverse atmospheric characteristics. Moreover, using the same Multi-Term LSM concept, the innovative multi-pixel concept has been implemented within the GRASP algorithm (Dubovik et al., 2011, 2021a). Under this concept, the inversion is performed simultaneously on a group of observations (e.g., satellite observations over different pixels in space and time). This allows for improved retrieval accuracy by applying additional a priori constraints on the spatial and temporal variability of retrieved parameters. This concept is vital for the implementation of the multi-platform retrieval because it enables synergetic retrieval of not fully coincident or not fully co-located observations.

The GRASP algorithm has already been successfully applied to observations from different spaceborne instruments. Extended aerosol characterization using the GRASP algorithm was demonstrated with PARASOL measurements (Popp et al., 2016; Chen et al., 2020; Schutgens et al., 2021). Application of GRASP to the Sentinel-3A/OLCI instrument showed AOD retrieval performance comparable to the MODIS dark target (DT) product (Chen et al., 2022b). Nevertheless, the studies also showed a strong dependence of GRASP/OLCI retrieval quality on a priori information about surface BRDF and reduced quality in the retrieval of extended properties like Angstrom Exponent (AE) and Single Scattering Albedo (SSA) (Chen et al., 2022b). Applying GRASP to Sentinel-5P/TROPOMI measurements demonstrated that the quality of extended aerosol and surface characterization can be significantly enhanced even from a single viewing instrument if conditions (ii), (iv), (v), and (vii) are fulfilled for the sensor and the retrieval algorithm (Litvinov et al., 2024; Chen et al., 2024a).

The GRASP multi-pixel retrieval strategy has already been used in various synergetic approaches. In particular, a synergetic approach was developed for sun-photometer and LIDAR ground-based measurements (Lopatin et al., 2013, 2021; Dubovik et al., 2021a). Synergetic retrieval from combined ground-based and satellite measurements was introduced and used to generate a surface reference database for validating satellite surface retrieval (GROSAT/GRASP approach; Litvinov et al., 2020, 2022, 2024). In Litvinov et al. (2021) and Chen et al. (2024b), a hybrid synergy with the GRASP algorithm takes advantage of the rich information content of TROPOMI measurements and the high spatial resolution of the PRISMA instrument (L2 to L1 hybrid synergy).

In this paper, we present a novel generalized synergetic approach implemented in the GRASP algorithm for retrieval of multi-mission spaceborne measurements – the SYREMIS/GRASP (SYnergetic REtrieval from Multi-MISsion instruments) approach – which can be robustly applied to present and future satellite observations. The concept was tested on two types of synergetic measurements: (i) synergy of Low Earth Orbiting (LEO) (polar-orbiting) satellites (LEO + LEO) and (ii) synergy of LEO and geostationary (GEO) satellites (LEO + GEO). The LEO + LEO synergy was implemented and tested on combined measurements from Sentinel-5P/TROPOMI, Sentinel-3A/OLCI, and Sentinel-3B/OLCI instruments (hereinafter also referred to as S5P/TROPOMI, S3A/OLCI, S3B/OLCI). The LEO + GEO synergy concept was applied to S5P/TROPOMI, S3A/OLCI, S3B/OLCI, and Himawari-8/AHI instruments.

In this paper, first, we describe the main principles of the developed SYREMIS/GRASP synergetic approach. Then, based on the validation results, the enhanced capabilities of the synergetic approach will be demonstrated, and its main drivers discussed.

2 SYREMIS/GRASP synergetic concept

The SYREMIS/GRASP synergetic approach was developed for currently operating polar-orbiting and geostationary satellites: S3A/OLCI, S3B/OLCI, S5P/TROPOMI, and Himawari-8/AHI. Figure 1 schematically demonstrates the general concept of combining multi-instrument spaceborne measurements within the SYREMIS/GRASP synergy. On the one hand, such multi-mission measurements allowed for testing the approach on actual aerosol events and evaluating the enhancement in aerosol characterization relative to already validated GRASP retrievals from S3A/OLCI, S5P/TROPOMI, and Himawari-8/AHI sensors alone (Chen et al., 2022b, 2024a; Litvinov et al., 2024; Li et al., 2025). On the other hand, it allowed filling gaps in detailed extended aerosol characterization, which have existed since the end of the POLDER-3/PARASOL mission in 2013 until the beginning of new polarimetric missions (HARP-2 and SPEX on board PACE, MetOp-SG 3MI, and others).

Figure 1 Schematic representation of the SYREMIS/GRASP multi-instrument measurements.

In general, such combined L1 measurements provide several advantages crucial for aerosol/surface characterization (Table 1): (i) better spectral coverage (from the UV to SWIR spectral range); (ii) improved temporal coverage with a few measurements per day in the LEO + LEO synergy and hourly (or better) temporal resolution in the LEO + GEO synergy; (iii) enhanced global coverage due to the combined measurements in the synergy; (iv) “quasi-multi-angular” measurements incorporating single-view observations from all instruments in Table 1, obtained at different illumination and observation geometries (solar and viewing zenith angles, azimuth angles). As a result, such combined synergetic measurements can be considered multi-angular or, strictly speaking, quasi-multi-angular, taking into account different measurement times for each observation angle.

Table 1 Multi-mission constellation for prototyped synergetic retrieval.

1*: https://user.eumetsat.int/resources/user-guides/sentinel-3-olci-level-1-data-guide (last access: 15 December 2025)
2*: https://sentinels.copernicus.eu/data-products (last access: 15 December 2025)
3*: https://www.eorc.jaxa.jp/ptree/userguide.html (last access: 15 December 2025)
4*: https://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/sample_netcdf.html (last access: 15 December 2025).

The advantages of combined L1 measurements are beneficial, provided that all different observations are simultaneous and/or co-located. However, this is usually not the case for observations from different satellite platforms. Therefore, the benefit of information complementarity in different observations is not straightforward. The application of the aforementioned multi-pixel concept enabled overcoming this issue. Specifically, a group of observations, including measurements from different satellites, is inverted simultaneously under a priori constraints on temporal and spatial variability of various atmospheric and surface reflectance parameters. As will be shown below, applying these constraints allows the propagation of information in the interpretation of all observations from different satellites.

Figure 2 illustrates the SYREMIS/GRASP multi-instrument Level 1 (L1) synergetic processing chain for three different instruments: S3A/OLCI, S3B/OLCI, and S5P/TROPOMI (marked with different colors). For each instrument, the cubic layers symbolize the measurements on a certain day, which represent the observed radiance at measurement geometry (solar and observation zenith angles, azimuth angles difference), at selected spectral bands, and a given spatial resolution. All measurements from different instruments are cloud-screened, regridded to a consistent spatial resolution, and then merged into a spatial–temporal multi-pixel block, which is used as L1 input for the SYREMIS/GRASP forward run and inversion. The synergetic SYREMIS product contains retrieved aerosol and surface properties for all spectral bands, derived at the time of the merged L1 measurements.

Figure 2 Illustration of SYREMIS multi-instrument processing chain.

Overall, the SYREMIS/GRASP synergetic concept is based on three main principles, which will be discussed in the next section:

• i.

  Harmonization and merging of L1 data from different instruments.

• ii.

  “Weighting” the used data sources (measurements and a priori constraints) according to their accuracies and information content.

• iii.

  Optimization of forward models and retrieval setup.

2.1 SYREMIS/GRASP synergetic measurements harmonization

All satellites from Table 1 differ significantly in their measurement capabilities, including spectral range and spatial coverage, spatial resolution, revisiting time. To exploit the extended information content of the synergy, measurements from each instrument should be properly harmonized and merged into a spatial–temporal multi-pixel group for simultaneous retrieval. In particular, this includes (i) spectral band selection within the synergy from different instruments; (ii) harmonization of spatial resolution and gridding; (iii) cloud masking merging.

Independent retrieval with the GRASP algorithm has already been applied to each instrument in Table 1, using 9 spectral bands for the OLCI sensor, 10 for TROPOMI, and 6 for the AHI instrument (Chen et al., 2022b, 2024a; Litvinov et al., 2024; Li et al., 2025; Table 2). These bands were also used in the SYREMIS/GRASP approach. In particular, 19 spectral measurements were included in the polar-orbiting LEO + LEO synergetic retrieval, and 24 bands in the LEO + GEO synergy (Table 2).

Table 2 Harmonized measurements from the SYREMIS synergy.

All instruments from the synergetic satellite constellation have different spatial resolutions (Table 1). To properly apply the multi-pixel constraints in the GRASP algorithm, all measurements were regridded to the same spatial grid. Following the approach used in GRASP/TROPOMI retrieval (Litvinov et al., 2024; Chen et al., 2024a), an equidistant cylindrical projection and WGS84 coordinate system with a spatial pixel resolution of 0.09° was used for all selected spectral bands, employing a bilinear regridding method (Table 2).

Since the measurements in the SYREMIS synergy are not collocated in time, cloud masking was applied independently for each instrument. In particular, cloud screening for TROPOMI in SYREMIS/GRASP processing is based on the S5P NPP-VIIRS cloud product, with approximately 500 m spatial resolution, as was done for GRASP/TROPOMI processing (Litvinov et al., 2024; Chen et al., 2024a). Similar to GRASP/OLCI retrieval (Chen et al., 2022b), IDEPIX cloud masking was applied to OLCI instruments (Wevers et al., 2022). Himawari-8/AHI cloud screening is based on Level 2 cloud products (Letu et al., 2018, 2020).

2.2 “Weighting” observations for the multi-instrument synergistic retrieval

The GRASP multi-pixel concept is a key methodological approach in realizing multi-platform inversion. It is implemented within the Multi-Term LSM methodology, as described by Dubovik et al. (2011, 2021a). Following this approach, the retrieval is conducted not for each selected observed satellite pixel, but rather for a set of observations collected over different locations at different time moments (spatial–temporal multi-pixel dataset block, Figs. 2 and 3).

Figure 3 Illustration of applying a multi-pixel approach for S3A/OLCI, S3B/OLCI, S5P/TROPOMI, and Himawari-8/AHI multi-platform retrieval.

Figure 3 illustrates an example of a spatial–temporal multi-pixel dataset block for 5×5 spatial pixels within a single selected day, filled with S3A/OLCI, S3B/OLCI, S5P/TROPOMI, and Himawari-8/AHI measurements. As one can see, the single-day dataset block consists of data from different sensors, performing observations at different times during the same day: S3A/OLCI and S3B/OLCI make observations in the morning and can overlap spatially for some pixels; S5P/TROPOMI is the afternoon satellite; Himawari-8/AHI is in geostationary orbit, providing multiple observations during a day over the same area.

The balanced measurement “weighting” in the spatial–temporal multi-pixel dataset block (Figs. 2 and 3) is crucial for synergetic retrieval. For example, as demonstrated in Litvinov et al. (2024) and Chen et al. (2024a), the S5P/TROPOMI instrument has the highest information content among all satellites from Table 1, allowing considerable enhancement of extended aerosol and surface characterization in comparison with OLCI and AHI instruments (Chen et al., 2022b; Li et al., 2025). This fact is accounted for by “weighting” measurements from different instruments in the SYREMIS/GRASP retrieval.

The GRASP algorithm is based on the statistical optimization concept (e.g., Dubovik et al., 2021a), where the contributions of different input data are balanced using known covariance matrices of each input dataset, including both satellite measurements and a priori data. Since those datasets are statistically independent (C=Iσ², where C is the covariance matrix, I is the diagonal unity matrix, and σ² is the variance (σ is the standard deviation) (Appendix A; Dubovik, 2004; Dubovik et al., 2021a)), the weights (importance of the measurements) are driven by known (or assumed) standard deviations in each dataset: the smaller the standard deviation of measurement fitting (the standard deviation of retrieval residual) required in the retrieval, the greater the “weight” of such measurements that can be assigned in the synergy. Within the framework of the Multi-Term LSM concept, the “weights” are quantitatively determined as ratios of ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ (variance of the first dataset) to ${\mathit{\sigma }}_i^{\mathrm{2}}$ (variance of the ith dataset). A detailed description of the data weighting Multi-Term LSM concept can be found in Dubovik et al. (2004, 2021a), and a brief discussion is also provided in Appendix A.

Table 3 shows the weighting of different instruments used in the multi-instrument synergy inversion. It can be seen that the S5P/TROPOMI has the highest weight among other observations. At the same time, it is important to emphasize that, in general, the weights are defined with respect to both measurements and the a priori datasets used. This aspect is discussed in the next section and Appendix A.

Table 3 Instrument “weighting” in SYREMIS LEO + LEO and LEO + GEO synergy setup: example over land.

2.3 Forward models and a priori constraints in the multi-instrument synergy

Application of the GRASP retrieval algorithm independently to OLCI, TROPOMI, and AHI single-view instruments, as well as to multi-angular polarimetric PARASOL measurements, showed good performance of the so-called aerosol “models” approach (Chen et al., 2020, 2022b, 2024a; Dubovik et al., 2021a; Hasekamp et al., 2024; Litvinov et al., 2024). In this approach, the total single-scattering characteristics of aerosols are represented as linear combinations of the characteristics of preselected aerosol components, i.e., assuming an external mixture of four different aerosol models: (1) fine absorbing, representing a climatological biomass burning aerosol model; (2) fine slightly absorbing, corresponding to climatological sulfate aerosols with introduced minor absorption; (3) coarse component representing mainly maritime aerosol type; and (4) coarse component for dust aerosol representation (Lopatin et al., 2021; Litvinov et al., 2024). The state vector for such an aerosol model in the GRASP inversion consists of the following characteristics: aerosol scale height, aerosol concentration, and aerosol model fractions (see, for example, the aerosol model parameters discussion in Litvinov et al., 2024).

The surface reflectance in the GRASP retrieval is usually described by the spectrally constrained renormalized Ross-Li BRDF model over land (Litvinov et al., 2010, 2011a, b, 2024) and a modified Cox and Munk model with accounting for water-leaving reflectance (Cox and Munk, 1954; Litvinov et al., 2024). These aerosol and surface reflectance forward approaches were found to be optimal also in the SYREMIS/GRASP synergetic retrieval for the satellite constellation from Table 1.

The strategic advantage of the GRASP inversion implemented via Multi-Term LSM statistically optimized fitting is that this concept allows the use of multiple a priori constraints in the retrieval. Such a priori constraints can be applied to all aerosol and surface characteristics and parameters (Dubovik et al., 2011, 2014, 2021a). In particular, the “single-pixel” a priori smoothness constraints can be used for limiting variability (and avoiding unrealistic oscillations) of, for example, aerosol column concentration, aerosol model fractions, spectral dependencies of surface BRDF parameters (the “in-pixel” smoothness constraints in Eq. A13). In addition to the “single-pixel” constraints, the “inter-pixel” constraints were used within the so-called “multi-pixel” retrieval when a large group of pixels was retrieved simultaneously and smoothness constraints (the “inter-pixels” smoothness constraints in Eq. A13) can be used for limiting temporal or spatial variability of parameters retrieved in different neighboring pixels (Dubovik et al., 2011, 2021a).

The successful application of the GRASP algorithm to different satellites demonstrated the crucial role of the single- and multi-pixels (or in-pixel and inter-pixel) smoothness constraints for stable and accurate characterization of both atmospheric aerosol and surface (Dubovik et al., 2011, 2014; 2021a; Chen et al., 2022b; Litvinov et al., 2024). A brief summary of the formal implementation of the “multi-pixel” retrieval is provided in Appendix A.

In comparison to the single-instrument retrieval, the SYREMIS LEO + LEO and LEO + GEO synergies have much better diurnal temporal resolution of measurements, when the temporal difference between S3A/OLCI, S3B/OLCI, S5P/TROPOMI, and especially between Himawari-8/AHI measurements can vary from several minutes to several hours within a single day. Moreover, the retrieval is performed on data accumulated in spatial–temporal blocks (Figs. 2 and 3) consisting of hundreds of satellite observations and covering about a month in the LEO + LEO and about two weeks in the LEO + GEO synergies (up to 150 and 200 temporal observations in the LEO + LEO and LEO + GEO synergies, respectively, Table 4). The temporal variability of aerosol and surface properties can be different within a few minutes, hours, or from day to day. Moreover, the temporal variability of aerosol properties is usually stronger than the surface ones. Therefore, the importance of using a priori temporal constraints for aerosol and surface parameters within a few minutes to several days and weeks significantly increases in the LEO + LEO and LEO + GEO configurations.

GRASP limits the temporal variability of the retrieved parameters by applying temporal a priori constraints, introduced as a priori estimates of the first derivatives with respect to time, approximated by finite differences (the differences between the parameters divided by the period over which they were observed; Dubovik et al., 2011, 2021a). The estimates are assumed to be normally distributed values with zero means and the standard deviations assigned on the basis of the expected temporal variability of retrieved parameters. The use of limitations applied to the finite differences is a very convenient way to ensure flexibility in temporal constraints. Indeed, such constraints allow much larger variability for parameters corresponding to more distant in time observations compared to the parameters from observations that are very close in time (full description can be found in Dubovik et al., 2011, 2021a).

However, in cases when observations are nearly simultaneous or very close in time, the finite differences may have very large values that can produce an imbalance in practical fitting. To avoid such difficulties, “temporal thresholds” on aerosol and surface characteristics are used in the SYREMIS/GRASP. The temporal threshold allows for avoiding very large values of the finite differences in cases when the combined synergistic observations are close to each other in time. Specifically, the very small time difference in the denominator of the finite difference is replaced by the specified “temporal threshold” value, allowing stronger temporal constraints on the temporal variability of the retrieved parameters within the threshold. This approach was in particular useful for constraining BRDF parameters, defined by the intrinsic properties of the surface (surface type, reflecting properties, topology), which are very stable in time, and, as a rule, do not change considerably during a day. The temporal thresholds applied to BRDF in the SYREMIS/GRASP approach allow using almost the same BRDF parameters for all synergetic measurements during the specified period. Similarly, the temporal thresholds are also useful for constraining the variability of aerosol model parameters (Table 4). Overall, in the SYREMIS/GRASP approach, different temporal thresholds were applied to different aerosol and surface properties, accounting for their different temporal dependence.

Table 4 Example of SYREMIS LEO + LEO and LEO + GEO smoothness constraints over land.

When the time interval between measurements in the synergy exceeds the specified temporal threshold, the multi-pixel temporal constraints (inter-pixel smoothness constraints on t-temporal variability in Eq. A13) in the SYREMIS/GRASP account for the actual temporal interval between the measurements, as was done in the single-instrument GRASP retrieval (Dubovik et al., 2011, 2021a; Chen et al., 2022b; Litvinov et al., 2024; Li et al., 2025). In addition, some spatial multi-pixel constraints were also used for aerosol, as described in Appendix A; however, these constraints played a rather minor role, as they were applied only within a 2-by-2-pixel spatial segment in the spatial–temporal multi-pixel block.

2.4 Remote sensing tests to optimize synergetic retrieval

A priori knowledge about single-instrument information content and measurement accuracy provides a good starting point for formulating the initial concept for the retrieval setups: instrument weighting, temporal, spatial, and spectral constraints in the synergy. Nevertheless, practice shows that optimizing and tuning this initial synergetic retrieval setup by conducting a series of remote sensing tests is usually desirable and necessary. In SYREMIS/GRASP, synergy optimization included several tests: (i) validation of retrieved aerosol properties over a limited number of AERONET stations; (ii) intercomparison of surface properties derived from synergetic ground-based and satellite measurements (based on the GROSAT/GRASP synergetic approach; Litvinov et al., 2020, 2022, 2024); and (iii) intercomparison of aerosol and surface synergetic retrieval with reference products from space-borne measurements over small regions ($\sim \mathrm{1000}\mathrm{km}\times \mathrm{1000}\mathrm{km}$).

Test (i) for the LEO + LEO and LEO + GEO synergies was first performed over a few AERONET stations (Holben et al., 1998; Giles et al., 2019) with different aerosol and surface properties (e.g., Mongu, Banizoumbou, Kanpur, Beijing) and then extended to 30 stations worldwide (e.g., Chen et al., 2022b) for the selected limited period (March–May 2019). The retrieval setups that showed the best AOD, AE, and SSA validation results were selected for further consideration. Test (ii) enabled intercomparison of surface BRDF parameters derived from different instruments using the GROSAT/GRASP synergetic approach, similar to Litvinov et al. (2022, 2024) and Chen et al. (2024a). Since these parameters are independent of illumination/observation geometries, they can be intercompared for different space-borne sensors, allowing identification of possible biases and their accounting for each instrument (Sect. 2.2). In Test (iii), MODIS and VIIRS (Levy et al., 2013; Sayer et al., 2018a; Schaaf and Wang, 2015a, b) aerosol/surface products were used for regional intercomparison with the SYREMIS/GRASP results.

The validation results against AERONET in Test (i) were evaluated using the following statistical characteristics: (i) Pearson correlation coefficient (R); (ii) Root Mean Square Error (RMSE); (iii) the number of data points (N) in the validation; and (iv) the number (and percentage) of data points satisfying the GCOS-based AOD criteria relative to AERONET. The GCOS-based criteria for AOD are based on the AOD requirements from Global Climate Observing System (GCOS-245, 2022) and their adaptation within the ESA aerosol-CCI project (Popp et al., 2016) for validation of AOD, derived from satellite measurements, against AERONET (de Leeuw et al., 2015): absolute difference in AOD is less than 0.04 or relative difference is less than 10 % (whichever is bigger).

The synergetic setup exhibiting the best validation statistical characteristics for aerosol and surface properties from Tests (i) and (ii), and which qualitatively agreed with the well-established aerosol/surface products (Test iii), was selected as a baseline approach for global synergetic retrieval.

The remote sensing tests, performed for the LEO + LEO and LEO + GEO synergies, showed that merging spectrally close measurements (for example, OLCI 442.5 nm and TROPOMI 440 nm, or OLCI 665 nm and TROPOMI 670 nm) does not provide the best results in Test (i). This can be due to different noise in the close spectral bands from different instruments, their different spectral response function, etc. Therefore, the final SYREMIS/GRASP synergetic data were prepared using 19 spectral measurements in the LEO + LEO synergy and 24 bands in the LEO + GEO synergy. The optimization test also showed the crucial role of balanced “weighting” of the multi-instrument measurements within the SYREMIS/GRASP synergetic approach. As discussed in Sect. 2.2 and Appendix A, instrument “weighting” in GRASP can be achieved by defining the standard deviation of the expected uncertainties, which determine weighting matrices in the GRASP multi-term LSM inversion scheme and, thus, the weights (or importance) of different data involved in the inversion. Experimenting with different values of the standard deviation, it was found that application of stronger requirements on the measurement fitting for the TROPOMI bands in the synergy resulted in much better validation results against AERONET (Test i) (Fig. 4). This can be explained by TROPOMI's higher information content (better spectral coverage and wider swath) in comparison with other instruments within the synergy. Specifically, whenever comparable values of the standard deviation were applied to all instruments and spectral bands, the validation against AERONET showed reduced performance relative to single-instrument retrieval (for example, in comparison with the GRASP/TROPOMI results; Litvinov et al., 2024; Chen et al., 2024a). Therefore, in the SYREMIS/GRASP synergetic approach, a few weighting groups with different requirements for the standard deviation of measurement fitting were associated with different instruments from the synergy (Table 3).

Figure 4 Comparison of the SYREMIS LEO + LEO AOD retrieval before (a) and after (b) harmonization, instrument “weighting”, and setting optimization. Each panel contains validation of AOD 550 nm vs. AERONET AOD 550 nm. In the scatter plot, the color of each datapoint indicates the two-dimensional probability density function (2D PDF) for AOD values (increasing PDF value corresponds to color gradient from blue to red). Below the scatter plot, the PDF AOD(SYREMIS) – AOD (AERONET) is provided for 4 AOD bins (see more details in Sect. 2.4). Left column of the panels: AOD product of all instruments in the synergy. Middle column of the panels: AOD product of S5p/TROPOMI extracted from the synergy. Right column of the panels: AOD product of S3A/OLCI and S3B/OLCI extracted from the synergy.

In addition, it was discovered that exchanging a few measurements between the weighting groups provided better retrieval results against AERONET. In particular, over land, the same standard deviation of 0.001 as for most TROPOMI channels (Table 3, group “a”) was assigned to OLCI bands 490 and 560 nm. Conversely, one TROPOMI band, 494 nm, was allocated to the OLCI bands' weighting group with a required standard deviation of 0.05 (Table 3, group “b”). Such band inter-exchange was found to be useful for improving OLCI retrieval in the LEO + LEO synergy without loss of accuracy in the SYREMIS/TROPOMI retrieval. This allows transferring information content from TROPOMI to OLCI and vice versa. The inter-exchange of selected bands between OLCI and AHI measurements also improved the consistency of retrieval for all instruments in the LEO + GEO synergy (Table 3). The optimal weighting of the different measurements used in the LEO + LEO and LEO + GEO SYREMIS synergies, according to the applied requirements on the standard deviation, is presented in Table 3.

As one can see, the “weighting” of the LEO + GEO is much more complicated than that of the LEO + LEO synergy. In particular, two additional groups of weights were added: group “c” with “0.01”, applied mainly to AHI instruments to account for diurnal aerosol variability from geostationary measurements; and group “d” with “0.05”, to reduce the “weight” of SWIR AHI measurements, that allowed improving the performance of the remote sensing tests.

The values of the “weights” in Table 3 are related to the specific inversion approach of the GRASP algorithm applied to the instruments listed in Table 1. Moreover, due to the ill-posed character of the inversion, the “weighting” from Table 3 is not unique even within the GRASP algorithm since other combinations of its values are possible with quite comparable retrieved results. The values provided in Table 3 demonstrate a general tendency for the considered synergies: the “weight” of S5P/TROPOMI measurements should be higher compared to other instruments.

Besides satellite measurements' “weighting”, GRASP “single-” and “multi-pixel” a priori constraints (Sect. 2.3, Appendix A) also play a key role in the successful realization of the SYREMIS/GRASP approach. Overall, the spectral and spatial constraints used in GRASP single-instrument retrieval (Chen et al., 2022b; Litvinov et al., 2024; Li et al., 2025) showed good performance in the optimization tests and were subsequently used in synergy.

Due to the merging of measurements from different sensors, the synergetic L1 input data blocks have much better temporal resolution compared to any single instrument from the synergy (Figs. 2 and 3). This affects the optimal temporal constraints applied to aerosol and surface parameters within the SYREMIS/GRASP approach. Properly selected temporal thresholds and smoothness constraints substantially increase the number of quasi-multi-angular measurements in the synergy, a factor crucial for BRDF parameter retrieval and the distinguishing of atmosphere and surface signals. The optimized temporal thresholds and smoothness constraints for aerosol and surface parameters in the SYREMIS/GRASP LEO + LEO and LEO + GEO synergies are presented in Table 4. As one can see, they account for known tendencies: surface BRDF temporal variability is generally smaller than the temporal variability of aerosol properties (stronger thresholds and constraints on the surface BRDF parameters than on the aerosol ones in Table 4); the temporal variability of the aerosol concentration is usually larger than that of the aerosol microphysics, which results in stronger constraints on the aerosol model fraction parameter in the SYREMIS/GRASP approach (Table 4).

The temporal thresholds applied to BRDF parameters in the SYREMIS approach considerably limit the variability of the BRDF parameters for all synergetic measurements during the specified period. They can change depending on pixel latitude due to different satellite overpass times at low and high latitudes. For the considered LEO + LEO synergetic satellite constellation, several hours (for example, ±6 h) are sufficient to account for a stable surface within the day globally.

LEO + GEO synergetic measurements cover a wider range of observation/illumination geometries (solar and viewing zenith and azimuth angles) than LEO + LEO. It has already been demonstrated that different BRDF models used in remote sensing perform well only for limited geometries (the limitations of different BRDF models in remote sensing are studied, for example, in Litvinov et al., 2011a, b). Due to this known limitation of the BRDF models, the optimized temporal thresholds for the LEO + GEO synergy are relaxed to a few hours, and the constraints on temporal variability rely mainly on temporal smoothness constraints in the GRASP algorithm, allowing the retrieved properties to change smoothly in time.

Similar to Table 3, due to the ill-posed nature of the inversion problem, the values provided in Table 4, strictly speaking, are not unique, since quite similar retrieval performance of aerosol and surface properties can be obtained within a certain range of the constraints. Nevertheless, the chosen parameters are generally within the optimal range in the sense that they adequately reflect the tendencies in temporal dependencies of aerosol and surface properties, adapted to the information content in the LEO + LEO and LEO + GEO synergies.

Figure 4 shows the SYREMIS LEO + LEO retrieval results obtained before (Fig. 4a) and after (Fig. 4b) harmonization, instrument “weighting”, “single-” and “multi-pixel” setup optimization. The results are intercompared with AERONET AOD at 550 nm. In this and other validation figures below for AOD, AE, and SSA, AERONET Level 2.0 datasets from the Version 3 Direct Sun and Inversion Algorithms were used (Holben et al., 1998; Dubovik and King, 2000; Dubovik et al., 2000; Giles et al., 2019; Sinyuk et al., 2020).

In Fig. 4 (as well as in other similar validation plots below), the color of each data point indicates the AOD Probability Density Function (PDF) value (an increasing PDF value corresponds to a color gradient from blue to red). The straight dotted line represents an ideal one-to-one correspondence in AOD. The straight red line is the linear fit to the data distribution, with the slope and intercept marked in the top-left corner legend. The grey envelope surrounding the one-to-one line in the scatter plot indicates the AOD GCOS-based envelope discussed above: an absolute difference in AOD is less than 0.04, or a relative difference is less than 10 % (whichever is bigger) (GCOS-245, 2022; de Leeuw et al., 2015). Below the scatter plots, the PDF (“Probability” in Fig. 4) of the AOD difference (AOD(SYREMIS/GRASP) – AOD(AERONET) in Fig. 4) is provided for four AOD intervals (bins). “Total” (black color in Fig. 4) corresponds to all AOD cases; “<0.2” (red): AOD<0.2; “[0.2,0.7]” (blue): AOD is between [0.2,0.7]; and “>0.7” (green in Fig. 4): AOD>0.7 cases. In Fig. 4, the bias (“Avg”), standard deviation (“StD”), Root Mean Square Error (“RMSE”), the total number of pixels in the validation (N), and the number and percentage of pixels satisfying GCOS-based requirements (“GCOS”) are indicated as a legend at the top of the PDF plots.

Table 5 lists the AOD validation statistical parameters for the SYREMIS LEO + LEO synergy extracted from Fig. 4. In particular, after optimizing the retrieval settings, “GCOS” (Total) improved from 34.2 %–47.9 %; bias (Total) decreased from 0.04 to −0.02; the correlation coefficient improved from 0.89–0.9; and the RMSE decreased from 0.16–0.14. The biggest improvements are for the low AOD cases (AOD<0.2), containing the majority of validation pixels: “GCOS” improved from 38.8 %–61 %, and “bias” decreased from 0.06–0.01. For moderate AOD ($\mathrm{0.2}<\text{AOD}<\mathrm{0.7}$), “GCOS” improved from 25.2 %–26.3 %, and “bias” changed from 0.04 to −0.04. For high AOD (>0.7) cases, “GCOS” remained flat around 35 %, and “bias” changed from −0.06 to −0.09. One can observe substantial improvement in the validation statistical characteristics, emphasized by the bold fonts in Table 5, relative to AERONET after properly accounting for the information content of the instruments, measurement accuracy, and adjustment of the retrieval approach to the synergetic multi-instrument constellation.

Table 5 Summary of the SYREMIS LEO + LEO AOD accuracy statistics for the two tests (Fig. 4a) “Before” and Fig. 4b “After”) in Fig. 4. The statistical parameters in the column headers are the same as in Fig. 4. “GCOS” and “Bias” are presented for each AOD bin (“Total”, “AOD<0.2”, “$\mathrm{0.2}\le \text{AOD}\le \mathrm{0.7}$”, “AOD>0.7”), same as in Fig. 4 histogram. The values in bold indicate better performance.

3 Validation and inter-comparison of the synergetic product

The validation of the SYREMIS/GRASP processing for LEO + LEO and LEO + GEO synergies was performed against AERONET and intercompared with VIIRS and MODIS aerosol and surface products (Schaaf et al., 2002; Schaaf and Wang, 2015; Hsu et al., 2003, 2019; Sayer et al., 2018a, b). The validation criteria are described in Sect. 2.4 for the optimization remote sensing test and were used for the GRASP/TROPOMI retrieval evaluation (Litvinov et al., 2024; Chen et al., 2024a).

3.1 SYREMIS/GRASP LEO + LEO synergy performance vs. AERONET

The validation results for the synergetic SYREMIS/GRASP LEO + LEO retrieval against the global AERONET stations for March, April, and May 2019 are presented in Figs. 5–7. Figures 5a–7a show the validation for all instruments in the synergy (SYREMIS LEO + LEO). Figures 5b–7b present the data associated with S5P/TROPOMI measurements (data extracted from the synergy at the time of TROPOMI measurements, indicated as “SYREMIS LEO + LEO: S5P/TROPOMI” in Figs. 5b–7b). Figures 5c–7c contain the data associated with S3/OLCI measurements (data extracted from the synergy at the time of S3A/OLCI and S3B/OLCI measurements, indicated as “SYREMIS LEO + LEO: S3/OLCI” in Figs. 5c–7c).

Figure 5 Validation of AOD 550 nm with global AERONET data over land and ocean for SYREMIS LEO + LEO synergy global processing in 2019 March, April, and May. Left column of panels: AOD product of all instruments in the synergy; middle column of panels: AOD product of TROPOMI extracted from the synergy; right column of panels: AOD product of S3A/OLCI and S3B/OLCI extracted from the synergy. The details of the scatter and PDF plots are explained in Fig. 4 and Sect. 2.4. The statistical metrics of each panel are summarized in Table 6.

Figure 6 The same as in Fig. 5, but for Angstrom Exponent validation vs. AERONET. Similar to AOD validation plots in Figs. 4 and 5, in the scatter plot, the color of each datapoint indicates the 2D probability density function value (increasing PDF value corresponds to color gradient from blue to red). In the legend in the top left corner of the scatter plots, the slope and intercept of the linear fit of the data points, the correlation coefficient R, and RMSE of the data points, and the number of data points N are indicated. Below each scatter plot, there is a corresponding histogram showing the error distribution, and in the legend over the top of the histogram, the bias between SYREMIS and AERONET AE (“Avg”), the standard deviation (“SD”), and the number of data points are indicated. The statistical metrics of each panel are summarized in Table 7.

Figure 7 The same as in Fig. 6 (AE) but for Single Scattering Albedo (SSA) 550 nm validation vs. AERONET. The statistical metrics of each panel are summarized in Table 7.

To remove the outliers from the retrieval and ensure consistency in the presented results, a filtering similar to the GRASP/TROPOMI quality assurance flag (Litvinov et al., 2024) was applied in Figs. 5–7. In particular, over land, the validation was done for pixels satisfying the following conditions: (i) the relative residual of fitting was less than 0.03 (3 %); (ii) the standard deviation (σ_AOD) of AOD (670 nm) within the 3×3 pixel window was less than 0.05, or the relative standard deviation σ_AOD (670 nm)$/$AOD (670 nm) was less than 0.15; (iii) the number of valid pixels within the 3×3 pixel window was greater than or equal to 5. Over the ocean, the filtering conditions were relaxed, taking into account the overall smaller total reflectance and global AOD values compared to pixels over land: (i) the relative residual was less than 0.1; (ii) σ_AOD (670 nm) < 0.05. For the Angstrom Exponent (AE) and SSA, an additional filter was applied when only pixels with AOD(550 nm) > 0.2 and AOD(550 nm) > 0.3, respectively, were used in the validation.

In general, the synergetic SYREMIS/GRASP retrieval shows good correspondence to AERONET, with a high percentage (∼54.3 % over land and ∼65 % over ocean) of fulfillment of GCOS-based (GCOS-245, 2022) requirements for AOD. The values of the statistical characteristics (Pearson correlation coefficient, Root Mean Square Error (RMSE), bias, etc.) indicate the high quality of the retrieval, as shown in Figs. 5–7.

The corresponding SYREMIS LEO + LEO AOD (550 nm) validation statistics are summarized in Table 6. As one can see, the SYREMIS all-instruments AOD (550 nm) retrievals (combining S5P/TROPOMI, S3A/OLCI, and S3B/OLCI, as presented in Fig. 5a and d) show good agreement with AERONET observations. Over land (Fig. 5a), the validation yields the correlation coefficient (R) of ∼0.89, RMSE of ∼0.094, and GCOS-based compliance fraction of 54.3 %. Over the ocean (Fig. 5d), the R value is ∼0.88, the RMSE is ∼0.060, and the GCOS-based compliance fraction is 65.0 %.

These AOD statistical characteristics are consistent with the validation results for all individual instruments extracted from the SYREMIS LEO + LEO retrievals, as shown in Fig. 5b and c (for S5P/TROPOMI and S3/OLCI over land) and Fig. 5e and f (for S5P/TROPOMI and S3/OLCI over ocean), respectively. Moreover, all retrieved extended characteristics (AE and SSA) are also consistent with each other across all instruments in the synergy. Specifically, the validation results for AE and SSA from the SYREMIS LEO + LEO retrievals are presented in Figs. 6 (AE) and 7 (SSA), respectively, and summarized in Table 7. AE from the SYREMIS LEO + LEO (Fig. 6a and d), as well as from the TROPOMI (Fig. 6b and e) and OLCI (Fig. 6c and f) data extraction, demonstrates similarly good agreement with AERONET, with correlation coefficients (R) ranging from 0.70–0.80 and RMSE values of 0.45–0.48 over land. Over the ocean, R values range from 0.60–0.75 with RMSEs of 0.40–0.50. For the SYREMIS LEO + LEO (Fig. 7a and d), as well as for TROPOMI (Fig. 7b and e) and OLCI (Fig. 7c and f) data extraction, the SSA (550 nm) validations against AERONET show that the RMSEs are approximately 0.03–0.035 over land and 0.05–0.07 over ocean, indicating, in general, good agreement taking into account the limited available matchups of SSA under moderate and high AOD conditions.

Table 6 Summary of AOD 550 nm validation statistics for the SYREMIS LEO + LEO global processing validation in Fig. 5. Each row in the table corresponds to one validation plot in Fig. 5.

Table 7 Summary of AE and SSA (550 nm) validation statistics for the SYREMIS LEO + LEO global processing validation in Figs. 6 and 7. Each row in the table corresponds to one validation plot in Figs. 6 and 7.

3.2 SYREMIS/GRASP LEO + LEO inter-comparison with GRASP single instrument retrieval over AERONET

Figures 8–13 demonstrate the added value of the synergetic approach, comparing validation results against AERONET of AOD, AE, and SSA obtained from the SYREMIS/GRASP LEO + LEO retrieval and from the corresponding GRASP single-instrument retrieval. Specifically, Figs. 8a–10a present the validation results from the SYREMIS LEO + LEO synergy associated with S5P/TROPOMI measurements (data extracted from the synergy at the time of TROPOMI measurements: “SYREMIS LEO + LEO: S5P/TROPOMI” in Figs. 8–10). Figures 8b–10b show the validation obtained from the GRASP/TROPOMI single-instrument retrievals (Litvinov et al., 2024; Chen et al., 2024a). Similarly, Figs. 10–13 present the validation for OLCI data extracted from the SYREMIS LEO + LEO synergy (“SYREMIS LEO + LEO: S3A/OLCI” in Figs. 11a–13a) and for the GRASP/OLCI-A single-instrument retrieval in Figs. 11b–13b (Chen et al., 2022). The corresponding intercomparison of validation statistics is summarized in Tables 8–11. In general, one can observe a higher number of retrievals passing the filtering criteria (mainly the residual filter criteria, Sect. 3.1) in SYREMIS LEO + LEO products compared to the GRASP single-instrument retrievals, indicating increased quality of the retrieval in the synergy. This improvement may be attributed to enhanced aerosol–surface signal decoupling in the synergy, resulting in smaller spectral fitting residuals. In addition, the GRASP/TROPOMI single-instrument AOD (550 nm) retrieval exhibits some scatter, with an RMSE of 0.139 and a GCOS compliance fraction of 48.4 % over land, as shown in Fig. 8b. These metrics are improved in the SYREMISLEO + LEO: S5P/TROPOMI retrievals, which achieve an RMSE of 0.090 and a GCOS fraction of 53.9 %, as shown in Fig. 8a. For the extended aerosol characteristics (AE and SSA), the SYREMIS LEO + LEO: S5P/TROPOMI retrievals show validation statistics that are generally comparable to those from the GRASP/TROPOMI single-instrument retrievals, as presented in Figs. 9 and 10, respectively. On the other hand, a clear improvement is observed for the SYREMIS LEO + LEO: S3A/OLCI aerosol data compared to the GRASP/OLCI-A single-instrument retrievals. Notably, the number of available AOD (550 nm) matchups with AERONET nearly doubles both over land and ocean (Fig. 11a and b). Additionally, the GCOS compliance fraction over land increases from 42.5 %–54.5 % (Fig. 11a and b), comparable to the SYREMIS LEO + LEO: S5P/TROPOMI. Although the validation statistics over the ocean show a slight decrease, as seen in Fig. 11, the number of successful inversions increases by a factor of three (162–543), and more high-AOD values pass the filtering criteria with 64.5 % high GCOS fulfillment, indicating enhanced sensitivity and broader retrieval capability in challenging conditions. For the extended aerosol characteristics (AE and SSA), as shown in Figs. 12 and 13, the validation results for the SYREMIS LEO + LEO: S3A/OLCI datasets are considerably improved in comparison with the single instrument GRASP/OLCI-A retrieval and reached a level comparable with those from the GRASP/TROPOMI single-instrument retrievals. These findings further support the conclusion that S5P/TROPOMI remains the dominant contributor to the information content in the SYREMIS LEO + LEO synergistic retrievals.

Figure 8 AOD validation from (a) SYREMIS LEO + LEO: TROPOMI and (b) GRASP/TROPOMI vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 8.

Figure 9 AE validation from (a) SYREMIS LEO + LEO: TROPOMI and (b) GRASP/TROPOMI vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 9.

Figure 10 SSA validation from (a) SYREMIS LEO + LEO: TROPOMI and (b) GRASP/TROPOMI vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 9.

Figure 11 AOD validation from (a) SYREMIS LEO + LEO: S3A/OLCI and (b) GRASP/OLCI-A vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 10.

Figure 12 AE validation from (a) SYREMIS LEO + LEO: S3A/OLCI and (b) GRASP/OLCI-A vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 11.

Figure 13 SSA validation from (a) SYREMIS LEO + LEO: S3A/OLCI and (b) GRASP/OLCI-A vs. AERONET over land and ocean in March, April, and May 2019. The statistical metrics of each panel are summarized in Table 11.

Summarizing intercomparison results, one can see better performance of the SYREMIS LEO + LEO: S5P/TROPOMI AOD data extract in comparison with the GRASP/TROPOMI, while the extended properties (AE and SSA) are of similar quality in both (synergy and single-instrument retrieval) cases. The biggest improvement is observed for the OLCI instrument, where the performance of all retrieved parameters in SYREMIS LEO + LEO: S3A/OLCI data extract is essentially improved compared to the single-instrument GRASP/OLCI-A retrieval. In particular, all retrieved characteristics from the synergy (AOD, AE, and SSA) are much better than those from the single instrument GRASP/OLCI-A retrieval, and the number of pixels passed through the quality filter considerably increased. In general, SSA from the SYREMIS LEO + LEO: S3/OLCI dataset is of the same quality as from the SYREMIS/LEO + LEO: S5P/TROPOMI, providing aligned retrieval for all instruments from the synergy.

The presented results clearly show that properly combining measurements from the S3A/OLCI, S3B/OLCI, and S5P/TROPOMI according to information content and accuracy, the synergetic SYREMIS/GRASP approach considerably improves the retrieval of the extended aerosol characterization in comparison to the single instrument retrieval. It enables transfer of information content from one instrument to another, which results in consistent aerosol characterization from non-coincident diverse satellite measurements. As a result, the quality of the retrieval for all instruments in the SYREMIS LEO + LEO synergy is comparable to or better than the quality from the GRASP/TROPOMI single-instrument retrieval and much better than from the GRASP/OLCI single-instrument retrieval. This can be explained by the fact that S5P/TROPOMI measurements have a wider spectral range and swath (consequently, better temporal resolution) than S3/OLCI. They are also known for their high radiometric accuracy (Ludewig et al., 2020; Tilstra et al., 2020). All these are crucial for atmosphere and surface signals differentiation and enhancement of the retrieval as demonstrated in Litvinov et al., 2024 and Chen et al., 2024a, and emphasized in the presented SYREMIS/GRASP synergetic results.

Table 8 Intercomparison of AOD (550 nm) AERONET validation statistics from SYREMIS LEO + LEO: S5P/TROPOMI and GRASP/TROPOMI single-instrument retrievals over land and ocean (Fig. 8). The best performing metric is indicated in bold.

Table 9 Intercomparison of AE and SSA (550 nm) AERONET validation statistics from the SYREMIS LEO + LEO: S5P/TROPOMI and GRASP/TROPOMI single-instrument retrievals over land and ocean (Figs. 9 and 10). The best performing metric is indicated in bold.

Table 10 Intercomparison of AOD (550 nm) AERONET validation statistics from SYREMIS LEO + LEO: S3A/OLCI and GRASP/OLCI single-instrument retrievals over land and ocean (Fig. 11). The best performing metric is indicated in bold.

NA: not available.

Table 11 Intercomparison of AE and SSA (550 nm) AERONET validation statistics from the SYREMIS LEO + LEO: S3A/OLCI and GRASP/OLCI single-instrument retrievals over land and ocean (Figs. 12 and 13). The best performing metric is indicated in bold.

3.3 SYREMIS/GRASP LEO + GEO synergy performance vs. AERONET

Validation of AOD, AE, and SSA from the SYREMIS LEO + GEO synergy over all AERONET stations within Himawari-8/AHI full scan is presented in Figs. 14–16 and summarized in Tables 12 and 13. Similar to the LEO + LEO synergy, the results are presented for all instruments involved in the SYREMIS/GRASP synergy, as well as for the data extracted for each specific instrument (SYREMIS LEO + GEO: S5P/TROPOMI, SYREMIS LEO + GEO: S3/OLCI (OLCI-A and OLCI-B), and SYREMIS LEO + GEO: Himawari-8/AHI). The retrievals of the AOD, AE, and SSA are highly comparable across all instruments in the synergy. For example, the AOD (550 nm) validation with AERONET shows the correlation coefficients around 0.92, the RMSE values near 0.15, and the GCOS-based fraction of approximately 40 %. While the RMSE and GCOS-based fraction are slightly worse than the global SYREMIS LEO + LEO AOD (550 nm) validation results presented in Fig. 5, this difference can be attributed to the Himawari-8/AHI full scan region, which includes a higher proportion of high-AOD cases and thus poses more challenging retrieval conditions. Moreover, for all instruments, AE (with RMSE ∼0.3) and SSA (with RMSE<0.05) are of similar quality to AE and SSA from GRASP/TROPOMI (Litvinov et al., 2024) and the SYREMIS LEO + LEO synergy (Figs. 6 and 7).

Figure 14 AOD 550 nm validation against AERONET data over land for SYREMIS LEO + GEO synergetic retrieval in 2019 March, April, and May. (a) Validation of AOD of all instruments in the synergy; (b) validation of AOD of S5p/TROPOMI extracted from the synergy; (c)  validation of AOD of S3/OLCI extracted from the synergy; (d) validation of AOD of Himawari-8/AHI extracted from the synergy. The statistical metrics are summarized in Table 12.

Figure 15 Angstrom Exponent (AE) validation against AERONET data over land for SYREMIS LEO + GEO synergetic retrieval in 2019 March, April, and May. (a) Validation of AE of all instruments in the synergy; (b) validation of AE of S5P/TROPOMI extracted from the synergy; (c) validation of AE of S3/OLCI extracted from the synergy; (d) validation of AE of Himawari-8/AHI extracted from the synergy. The statistical metrics are summarized in Table 13.

Figure 16 Single Scattering Albedo (SSA) 550 nm validation against AERONET data over land for SYREMIS LEO + GEO synergetic retrieval in 2019 March, April, and May. (a) Validation of SSA of all instruments in the synergy; (b) validation of SSA of S5p/TROPOMI extracted from the synergy; (c) validation of SSA of S3/OLCI extracted from the synergy; (d) validation of SSA of Himawari-8/AHI extracted from the synergy. The statistical metrics are summarized in Table 13.

Table 12 Summary of AOD validation statistics for the SYREMIS LEO + GEO processing in Fig. 14. Each row in the table corresponds to one validation plot in Fig. 14.

Table 13 Summary of Angstrom Exponent and Single Scattering Albedo validation statistics for the SYREMIS LEO + GEO processing in Figs. 15 and 16. Each row in the table corresponds to one validation plot in Figs. 15 and 16.

In general, the validation and inter-comparison results show that, similar to the LEO + LEO synergy, richer information content from the S5P/TROPOMI propagates to other instruments of the SYREMIS LEO + GEO synergy (S3A/OLCI, S3B/OLCI, and Himawari-8/AHI), improving extended aerosol characteristics such as AE and spectral SSA. At the same time, AOD is improved for all instruments in the SYREMIS LEO + GEO synergy, including S5P/TROPOMI, due to the additional spatial, temporal, and spectral information content in the synergetic measurements. Moreover, one of the crucial advantages of the LEO + GEO synergy is the high temporal resolution of the extended aerosol characterization (AOD, AE, SSA, etc.). In particular, the considered SYREMIS LEO + GEO synergy provides diurnal variability of aerosol with ∼1 h temporal resolution, which allows for monitoring aerosol transport, air quality, and can be used in atmospheric dynamics studies. The observed unique advantages in capturing the diurnal variability of aerosol optical-microphysical properties from the LEO + GEO SYREMIS product will be reported and discussed in a separate study.

3.4 SYREMIS/GRASP aerosol and surface products global intercomparison

In this section, we present a preliminary intercomparison of aerosol and surface products, obtained from the SYREMIS/GRASP synergetic retrieval, with GRASP/TROPOMI, VIIRS, and MODIS aerosol and surface products (Hsu et al., 2013; 2019; Sayer et al., 2018a; Schaaf and Wang, 2015a, b; Chen et al., 2022b; Litvinov et al., 2024). A more detailed analysis of the global SYREMIS product will be performed in separate studies.

Figure 17 Global map of monthly mean AOD at 550 nm of SYREMIS/TROPOMI, GRASP/TROPOMI, and VIIRS/DB in March 2019. The difference between the three products is also presented (“VIIRS/DB-TROPOMI/GRASP”, “VIIRS/DB-SYREMIS/TROPOMI”, “TROPOMI/GRASP-TROPOMI/SYREMIS”).

Figure 17 shows the global maps of the March 2019 monthly mean AOD at 550 nm for the SYREMIS LEO + LEO: S5P/TROPOMI (TROPOMI data extraction from the synergy, indicated as SYREMIS/TROPOMI in Fig. 17), GRASP/TROPOMI (single-instrument GRASP retrieval from TROPOMI), VIIRS/DB, and their differences. The SNPP/VIIRS AERDB_L2_VIIRS_SNPP aerosol product (Hsu et al., 2013, 2019) was used in the intercomparison. All products were regridded to the same 0.2° global equi-rectangular grid for pixel-to-pixel intercomparison.

Globally, SYREMIS/TROPOMI AOD values show good qualitative agreement with single-instrument GRASP/TROPOMI and VIIRS products. Nevertheless, the quantitative intercomparison (see the maps of the differences “VIIRS/DB-SYREMIS/TROPOMI” and “TROPOMI/GRASP-TROPOMI/SYREMIS” in Fig. 17) shows that the synergetic SYREMIS/TROPOMI AOD is lower than the GRASP/TROPOMI and VIIRS products over: (i) bright surfaces (e.g., over the Sahara and the Arabian Peninsula); (ii) the North Indian Subcontinent and the East coast of China, where AOD is quite large for all products; and (iii) the South polar regions, where AOD is very small for all products. In general, AOD products from the single instruments (GRASP/TROPOMI and VIIRS/DB) correspond to each other better than they correspond to the SYREMIS/GRASP AOD (see maps of the differences: “VIIRS/DB-TROPOMI/GRASP”, “VIIRS/DB-SYREMIS/TROPOMI”, and “TROPOMI/GRASP – TROPOMI/SYREMIS” in Fig. 17). A more detailed analysis of the global aerosol characterization with the synergetic approach will be the subject of separate studies.

Figure 18 Global intercomparison of SYREMIS/TROPOMI and MODIS Ross-Li BRDF 1st model parameter (isotropic) at “blue” band (440 and 470 nm) (a) and “red” band (670 and 660 nm) (b). The difference between the two products is presented as “MODIS-TROPOMI” at “blue” and “red” bands.

Figure 18 shows a qualitatively similar retrieval of the first (isotropic) BRDF parameter from the SYREMIS/TROPOMI data extraction compared to the MODIS product (the MODIS MCD43C1 Terra + Aqua BRDF Model Parameters Daily L3 Global 0.05Deg CMG product suite (Schaaf and Wang, 2015b) was used), regridded to the same 0.2° global equirectangular grid for pixel-to-pixel intercomparison. Depending on the region, SYREMIS/GRASP can provide a brighter or darker surface compared to MODIS. In the blue spectral band, SYREMIS/TROPOMI BRDF1 (the isotropic Ross-Li BRDF parameter) values are generally lower than MODIS BRDF1 values in the Saharan desert, the Middle East, and Tibet, and higher than MODIS BRDF1 in other regions. In the red spectral band, the difference between SYREMIS/TROPOMI and MODIS BRDF1 shows an opposite spatial distribution. In Central Asia and Tibet, MODIS BRDF1 is higher than SYREMIS/TROPOMI BRDF1 in both the 'blue' and 'red' spectral bands. SYREMIS and MODIS BRDF maps differ in spatial completeness over the Amazon and high-latitude regions due to differences in their cloud/snow masks.

Figure 19 Global intercomparison of SYREMIS/TROPOMI and MODIS Ross-Li BRDF 2nd (volumetric) (a) and 3rd (geometric) (b) model parameters.

Figure 20 Correlation of global BRDF second and third parameters between SYREMIS/TROPOMI and MODIS in March 2019. N is the number of pixels; R is the Pearson correlation coefficient.

The largest difference in the surface retrieval emerging from the synergetic approach can be found in the second (volumetric) and third (geometric) Ross-Li BRDF parameters, which describe the angular profile of the surface reflectance. As can be seen from Fig. 19, the values of these parameters derived by SYREMIS/GRASP over bright surfaces (e.g., over the Sahara) are much higher than those from the MODIS BRDF product. This can be clearly seen from the pixel-to-pixel correlation in Fig. 20. The stronger variability of the second and third BRDF parameters can be explained by the quasi-multi-angular measurements in the synergetic retrieval. In particular, the SYREMIS LEO + LEO accounts for S5P/TROPOMI, S3A/OLCI, and S3B/OLCI measurements within about one month, using temporal thresholds and “multi-pixel” temporal constraints of the GRASP algorithm (Table 4). This allowed for the accumulation of up to hundreds of synergetic measurements obtained at different observation/illumination geometries (zenith and azimuthal angles of the sun and satellites), which are accounted for in the GRASP multi-pixel LSM method (see Sect. 2). As a result, the synergetic retrieval of the angular features of BRDF (provided by the volumetric and geometric Ross-Li BRDF parameters) is expected to be much more comprehensive than any single-instrument retrieval, including a multi-day MODIS Terra+Aqua combination. Similar to the global aerosol product, synergetic surface BRDF retrieval on a global scale will be the subject of subsequent studies.

4 Conclusions

In this paper, we presented an approach for SYnergetic REtrieval from Multi-MISsion instruments (SYREMIS), which is based on the GRASP algorithm and was developed for enhanced retrieval of aerosol and surface properties. The approach was implemented and tested on (i) synergy from polar-orbiting satellites (the LEO + LEO synergy of the S5P/TROPOMI, S3A/OLCI, and S3B/OLCI instruments) and (ii) synergy of polar-orbiting and geostationary satellites (the LEO + GEO synergy of the S5P/TROPOMI, S3A/OLCI, S3B/OLCI, and Himawari-8/AHI instruments). The synergetic concept of the SYREMIS/GRASP approach is based on three main principles discussed in this paper (e.g., Sect. 2): (i) harmonization of multi-instrument L1 measurements, (ii) “weighting” the multi-instrument measurements, and (iii) optimization of the forward models and the retrieval setups.

Overall, it was demonstrated that by harmonizing multi-instrument measurements and properly balancing the “weights” of measurements from different sensors, the SYREMIS/GRASP approach allows for information transfer between all instruments in synergy. In combination with adjusted forward models and retrieval setups (spectral, spatial, and temporal constraints on the aerosol/surface variability), this results in increased performance of AOD, aerosol size, and absorption properties retrieval and more consistent surface BRDF characterization.

Performed remote sensing tests (e.g., Sect. 2.4) showed that the best SYREMIS/GRASP retrieval can be achieved when the “weight” of TROPOMI measurements in both LEO + LEO and LEO + GEO synergies is higher than the “weight” of OLCI and AHI instruments. TROPOMI is known for its rich information content in terms of spectral measurements (from the UV to SWIR spectral range) and wide swath (∼2600 km). In particular, UV measurements provide sensitivity to absorption properties; extended spectral coverage and a wide swath are crucial for aerosol size and type characterization, as well as for aerosol and surface signal differentiation (Litvinov et al., 2024). As shown in Sect. 3, within the considered synergetic satellite constellation, the TROPOMI instrument, providing the most information about aerosol and surface, serves as a “driver” of the SYREMIS/GRASP retrieval. Nevertheless, it was also demonstrated that the contribution of the other satellites in the SYREMIS LEO + LEO and LEO + GEO synergy is not negligible since they considerably extend the spectral, temporal, and spatial coverage of the TROPOMI instrument (e.g., Sect. 3). This allows us to conclude that the SYREMIS/GRASP approach, applied to synergetic multi-spectral, quasi-multi-angular measurements from multi-mission instruments, results in enhanced aerosol and surface BRDF characterization.

In particular, for all instruments from the LEO + LEO and LEO + GEO synergies, AOD performance against AERONET was increased. AE and SSA from the synergies were found to be comparable with the TROPOMI/GRASP single instrument retrieval and better than or of the same quality as previous OLCI/GRASP and AHI/GRASP retrievals (Sect. 3), supporting the conclusion that information about aerosol size and absorption/scattering comes mainly from the TROPOMI measurements.

SYREMIS/GRASP AOD inter-comparison with VIIRS and TROPOMI/GRASP shows similar global features, though observed regional differences require further analysis. Global inter-comparison of the retrieved surface properties showed good qualitative consistency of the SYREMIS/GRASP first BRDF parameter (isotropic Ross-Li parameter) with the MODIS parameter. Nevertheless, depending on the region, a brighter or darker surface may be retrieved. The SYREMIS/GRASP also shows differences in the volumetric and geometric Ross-Li BRDF parameters compared to the MODIS surface product.

This can be explained by the fact that SYREMIS multi-instrument quasi-multi-angular measurements provide more information about the angular dependence of surface reflectance compared to a single instrument with single-angle observations, as discussed in Sect. 3.4.

With properly applied temporal constraints, described in Sect. 2, the SYREMIS/GRASP retrieval allows for the derivation of consistent temporal variation in aerosol properties from all instruments within the synergy. In particular, the LEO + LEO synergy can provide extended aerosol properties several times per day, whereas the LEO + GEO synergy improves temporal resolution to one hour or better. Such extended aerosol characterization with high temporal resolution is required in air quality studies, for monitoring aerosol transport, aerosol-cloud interaction, etc.

The developed SYREMIS/GRASP synergetic approach is based on the fundamental principles of the multi-term LSM of the GRASP algorithm (e.g., Sect. 2 and Appendix A). Moreover, the formulated three main principles of the synergetic concept (e.g., Sect. 2) are quite universal and can be extended to future spaceborne missions, including synergy with multi-angular, multi-spectral, polarimetric measurements. In such advanced synergy with spaceborne polarimeters, more complex aerosol models can be used. This is expected to allow further enhancement of the retrieval of aerosol microphysical and chemical properties, as well as their temporal resolution.

Appendix A

The simultaneous solution for the selected data set corresponds to the solution of a joint system of equations that schematically can be represented as:

$$\begin{array}{}\text{(A1)}& \left\{\begin{array}{ll}{\mathit{f}}_{\mathrm{1}}^{\ast }& {\mathit{f}}_{\mathrm{1}}\left({\mathit{a}}_{\mathrm{1}}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathrm{1}}^{\ast }\\ {\mathit{f}}_{\mathrm{2}}^{\ast }=& {\mathit{f}}_{\mathrm{2}}\left({\mathit{a}}_{\mathrm{2}}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathrm{2}}^{\ast }\\ {\mathit{f}}_{\mathrm{3}}^{\ast }& {\mathit{f}}_{\mathrm{3}}\left({\mathit{a}}_{\mathrm{3}}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathrm{3}}^{\ast }\end{array}\right.=\mathit{f}\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{f}}^{\ast }.\end{array}$$

Here, the equations are written only for 3 pixels to illustrate the concept that can be extended to N_k pixels. ${\mathit{f}}_k^{\ast }$ corresponds to measurements associated with pixel k, f_k measurements simulation performed for the state vector a_k. Evidently, solving the system of independent equations (Eq. A1) simultaneously or consequently does not make any difference since each equation corresponds to a different pixel that depends on independent parameters. This is especially evident if the function ${\mathit{f}}_j^{\ast }={\mathbf{K}}_j{\mathit{a}}_i$:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\left\{\begin{array}{c}{\mathit{f}}_{\mathrm{1}}^{\ast }\\ {\mathit{f}}_{\mathrm{2}}^{\ast }\\ {\mathit{f}}_{\mathrm{3}}^{\ast }\end{array}=\left(\begin{array}{ccc}{\mathbf{K}}_{\mathrm{1}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_{\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& \mathrm{0}& {\mathbf{K}}_{\mathrm{3}}\end{array}\right)\left(\begin{array}{c}{\mathit{a}}_{\mathrm{1}}\\ {\mathit{a}}_{\mathrm{2}}\\ {\mathit{a}}_{\mathrm{3}}\end{array}\right)+\left(\begin{array}{c}\mathrm{\Delta }{\mathit{f}}_{\mathrm{1}}^{\ast }\\ \mathrm{\Delta }{\mathit{f}}_{\mathrm{2}}^{\ast }\\ \mathrm{\Delta }{\mathit{f}}_{\mathrm{3}}^{\ast }\end{array}\right)\right.\\ {\displaystyle }& {\displaystyle }=\mathit{f}\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{f}}^{\ast }\end{array}\end{array}$$

However, this system can be complemented by the system representing a priori constraints that provides some dependence between parameters a₁, a₂ and a₃:

$$\begin{array}{}\text{(A3)}& \left\{\begin{array}{ll}{\mathit{g}}_t^{\ast }& {\mathit{g}}_t({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}})+\mathrm{\Delta }{\mathit{g}}_t={\mathit{g}}_t\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{g}}_t\\ {\mathit{g}}_x^{\ast }=& {\mathit{g}}_x({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}})+\mathrm{\Delta }{\mathit{g}}_x={\mathit{g}}_x\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{g}}_x\\ {\mathit{g}}_y^{\ast }& {\mathit{g}}_y({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}})+\mathrm{\Delta }{\mathit{g}}_y={\mathit{g}}_y\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{g}}_y\end{array}\right.,\end{array}$$

where the a priori functional constraints g_t, g_x and g_y apply temporal and spatial dependencies between ${\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}}$.

Solving two systems (Eqs. A2 and A3) simultaneously changes the solution compared to the consequent solution of Eq. (A2). Since both measurements ${\mathit{f}}_{\mathrm{\dots }}^{\ast }$ and a priori ${\mathit{g}}_{\mathrm{\dots }}^{\ast }$ data include random errors Δ_…, the joint system (Eqs. A2 and A3) can be solved using Multi-Term LSM (Dubovik et al., 2021a) and corresponds to the minimum of the following quadratic form:

$$\begin{array}{}\text{(A4)}& \begin{array}{l}\mathrm{\Psi }\left(\mathit{a}\right)\\ ={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}\mathrm{\Delta }{\mathit{f}}_{\mathit{i}}^{\mathbf{T}}{\mathbf{C}}_i^{-\mathrm{1}}\mathrm{\Delta }{\mathit{f}}_i+{\sum }_{k=t,x,y}\mathrm{\Delta }{\mathit{g}}_{\mathit{k}}^{\mathbf{T}}{\mathbf{C}}_k^{-\mathrm{1}}\mathrm{\Delta }{\mathit{g}}_{\mathit{k}}\\ ={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}{\mathrm{\Psi }}_i\left({\mathit{a}}_i\right)+{\sum }_{k=t,x,y}{\mathrm{\Psi }}_k({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}})\\ ={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}{\mathrm{\Psi }}_i\left({\mathit{a}}_i\right)+{\sum }_{k=t,x,y}{\mathrm{\Psi }}_k\left(\mathit{a}\right)\\ ={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}{\mathrm{\Psi }}_i\left({\mathit{a}}_i\right)+{\mathrm{\Psi }}_{\text{inter-pixel}}\sim min\end{array}\end{array}$$

where $\mathrm{\Delta }{\mathit{f}}_{\mathit{i}}={\mathit{f}}_i^{\ast }-{\mathit{f}}_{\mathit{i}}\left({\mathit{a}}_i\right)$, C is the covariance matrix, $\mathbf{C}=\mathbf{I}{\mathit{\sigma }}^{\mathrm{2}},\mathbf{I}$ is the diagonal unity matrix, σ² is the variance, and σ is the standard deviation of the dataset.

Since the solution is related to the minimum of a quadratic form and not related to the absolute value of this minimum, Eq. (A4) can be rewritten using weighting matrices:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\displaystyle }\mathrm{\Psi }\left(\mathit{a}\right)& {\displaystyle }={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}{\mathit{\gamma }}_i\mathrm{\Delta }{\mathit{f}}_{\mathit{i}}^{\mathbf{T}}{\mathbf{W}}_i^{-\mathrm{1}}\mathrm{\Delta }{\mathit{f}}_{\mathit{i}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{k=t,x,y}{\mathit{\gamma }}_k\mathrm{\Delta }{\mathit{g}}_{\mathit{k}}^{\mathbf{T}}{\mathbf{W}}_k^{-\mathrm{1}}\mathrm{\Delta }{\mathit{g}}_{\mathit{k}}\sim min,\end{array}\end{array}$$

where W_… represent weighting matrices and γ_… represent so-called Lagrange multipliers, defined as

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{\displaystyle }{\mathbf{W}}_{\mathit{i}}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\{}{\mathbf{W}}_{\mathit{i}}{\mathit{\}}}_{\mathrm{11}}}}{\mathbf{W}}_{\mathit{f}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_i^{\mathrm{2}}}}{\mathbf{W}}_{\mathit{f}}\text{and}\\ {\displaystyle }{\mathit{\gamma }}_i=& {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}}\text{and}{\mathit{\gamma }}_k={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_k^{\mathrm{2}}}}.\end{array}\end{array}$$

The introduction of weighting matrices and Lagrange multipliers is a rather conventional approach (see discussion by Dubovik, 2004, Dubovik et al., 2021a) that is quite convenient since they allow explicitly expressing the weights (“importance”) of different data involved in the inversion.

Thus, using a priori constraints ${\mathit{g}}_k^{\ast }$ in Eq. (A3) allows for improving the accuracy of retrieval of multi-platform inversion. Formally, ${\mathit{g}}_k^{\ast }$ are usually defined by assuming the derivatives of temporal or spatial dependencies of retrieved parameters to zeros, i.e.

$$\begin{array}{rl}{\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{{\partial }^p{a}_n\left(x\right)}{\partial x^p}}\right)}^{\ast }\approx {\mathrm{\Delta }}_x^{\ast }={\mathrm{0}}_x^{\ast }+{\mathrm{\Delta }}_x,\\ {\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{{\partial }^p{a}_n\left(y\right)}{\partial x^p}}\right)}^{\ast }\approx {\mathrm{\Delta }}_y^{\ast }={\mathrm{0}}_y^{\ast }+{\mathrm{\Delta }}_y\text{and}\\ {\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{{\partial }^p{a}_n\left(t\right)}{\partial x^p}}\right)}^{\ast }\approx {\mathrm{\Delta }}_t^{\ast }={\mathrm{0}}_t^{\ast }+{\mathrm{\Delta }}_t\end{array}$$

(e.g. Dubovik et al., 2011, 2021a). Correspondingly, the multi-pixel constraints in Eq. (A3) are ${\mathit{g}}_k^{\ast }={\mathbf{0}}_k^{\ast }={\mathbf{S}}_k\mathit{a}+{\mathrm{\Delta }}_k$, where S_k is the numerical equivalent of corresponding kth derivatives, and it can be written:

$$\begin{array}{}\text{(A7)}& \begin{array}{l}{\mathrm{\Psi }}_{\text{inter-pixel}}={\sum }_{k=t,x,y}{\mathit{a}}^{\mathbf{T}}{\mathit{\gamma }}_k{\mathbf{\Omega }}_k\mathit{a}\\ \text{where}{\mathbf{\Omega }}_k={\mathbf{S}}_{\mathit{k}}^{\mathbf{T}}{\mathbf{W}}_k^{-\mathrm{1}}{\mathbf{S}}_k.\end{array}\end{array}$$

As shown before, the vector of retrieved parameters a includes the vectors of parameters retrieved for each pixel: ${\mathit{a}}_{\mathrm{1}}=({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},{\mathit{a}}_{\mathrm{3}})$ and each vector a_n includes the following retrieved characteristics of aerosol and surface:

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{\displaystyle }{\mathit{a}}_i^T=& {\displaystyle }({\mathit{a}}_{\mathrm{1}},\mathrm{\dots },{\mathit{a}}_{\mathrm{6}}{)}^{\mathrm{T}}\\ {\displaystyle }=& {\displaystyle }({\mathit{a}}_{C_{\mathrm{v}}},{\mathit{a}}_{\text{comp}},{\mathit{a}}_{\mathrm{h}},{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}{)}^{\mathrm{T}},\end{array}\end{array}$$

where:

$$\begin{array}{}\text{(A9)}& {\mathit{a}}_{\text{comp}}^{\mathrm{T}}={\left({\displaystyle \frac{{\mathit{a}}_{c_{\mathrm{v}\mathrm{1}}}}{{\mathit{a}}_{C_{\mathrm{v}}}}},{\displaystyle \frac{{\mathit{a}}_{c_{\mathrm{v}\mathrm{2}}}}{{\mathit{a}}_{C_{\mathrm{v}}}}},\mathrm{\dots },{\displaystyle \frac{{\mathit{a}}_{c_{\mathrm{v}n}}}{{\mathit{a}}_{C_{\mathrm{v}}}}}\right)}^{\mathrm{T}}\end{array}$$

where a_comp represents the unknown vectors of aerosol compositions, ${\mathit{a}}_{C_{\mathrm{v}}}$ and a_h represent the unknown vectors corresponding to total aerosol volume concentration and aerosol scale height. Over land, a_brdf,1, a_brdf,2 and a_brdf,3 represent the unknown vectors of the first, second, and third spectrally dependent surface Ross-Li BRDF parameters. Over ocean, a_brdf,1, a_brdf,2 and a_brdf,3 represent the unknown vectors of the ocean isotropic albedo, the fraction of Fresnel reflection, and the mean square facet slope as described previously by Litvinov et al. (2024). The inter-pixel constraints are usually different for different retrieved characteristics and, therefore, “inter-pixel spatial and temporal smoothness matrices” γ_kΩ_k have the following diagonal structure:

$$\begin{array}{}\text{(A10)}& {\mathit{\gamma }}_k{\mathbf{\Omega }}_k=\left(\begin{array}{cccc}{\mathit{\gamma }}_{k,\mathrm{1}}{\mathbf{\Omega }}_{k,\mathrm{1}}& \mathbf{0}& \mathrm{\dots }& \mathbf{0}\\ \mathbf{0}& {\mathit{\gamma }}_{k,\mathrm{2}}{\mathbf{\Omega }}_{k,\mathrm{2}}& \mathrm{\dots }& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathrm{\dots }& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathrm{\dots }& {\mathit{\gamma }}_{k,n}{\mathbf{\Omega }}_{k,n}\end{array}\right)\end{array}$$

where Lagrange multiplier γ_k,n and smoothness matrices Ω_k define inter-pixel constraints for each of retrieved characteristics: ${\mathit{a}}_{C_{\mathrm{v}}},{\mathit{a}}_{\text{comp}},a_{\mathrm{h}},{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}$.

It should also be noted that the a priori in retrieval is not only about a priori known inter-pixel relationships of the retrieved characteristics but also about known a priori knowledge variability of parameters within the pixels. Moreover, in each pixel we may have observations from different satellites, for example, S3A/OLCI, S3B/OLCI, S5P/TROPOMI, HIMAWARI-8/AHI, with different spectral observations. Specifically, the system of equations corresponding to observations for each ith pixel ${\mathit{f}}_i^{\ast }=\mathit{f}\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_i$ in Eq. (A2), consists of several lines:

$$\begin{array}{}\text{(A11)}& \left\{\begin{array}{l}{\mathit{f}}_{i,\mathrm{1}}^{\ast }={\mathit{f}}_{i,\mathrm{1}}\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,\mathrm{1}}\\ {\mathit{f}}_{i,\mathrm{2}}^{\ast }={\mathit{f}}_{i,\mathrm{2}}\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,\mathrm{2}}\\ \mathrm{\dots }\\ {\mathit{f}}_{i,N}^{\ast }={\mathit{f}}_{i,N}\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,N}\\ {\mathit{f}}_{i,\mathrm{1}}^a={\mathit{f}}_{i,\mathrm{1}}^a\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,\mathrm{1}}\\ {\mathit{f}}_{i,\mathrm{2}}^a={\mathit{f}}_{i,\mathrm{2}}^a\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,\mathrm{2}}\\ \mathrm{\dots }\\ {\mathit{f}}_{i,M}^a={\mathit{f}}_{i,M}^a\left({\mathit{a}}_{\mathit{i}}\right)+\mathrm{\Delta }{\mathit{f}}_{i,M}\end{array}\right..\end{array}$$

where ${\mathit{f}}_{i,n}^{\ast }$ represent measurements from different instruments or from the same instrument but with different accuracy, and ${\mathit{f}}_{i,m}^a$ represent a priori data that are usually defined as smoothness constraints defined as limitations on k-the derivatives, as ${\mathbf{0}}_k^{\ast }={\mathbf{S}}_k\mathit{a}+{\mathrm{\Delta }}_k$, or as direct a priori estimates ${\mathit{a}}_m^{\ast }={\mathit{a}}_m+{\mathrm{\Delta }}_{a,m}$, as suggested by Dubovik et al. (2011). Usually, in pixel a priori constraints are defined in the same way and with the same strength for all pixels. Therefore, the quadratic form Ψ_i,pixel(a_i) corresponding to ith pixel in Eq. (A3) can be written as follows:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\displaystyle }{\mathrm{\Psi }}_{i,}& {\displaystyle }={\sum }_{n=\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}}{\mathit{\gamma }}_{i,n}\mathrm{\Delta }{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}^{\mathbf{T}}{\mathbf{W}}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}^{-\mathrm{1}}\mathrm{\Delta }{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{m=\mathrm{1},\mathrm{\dots },M_{\mathrm{\Delta }}}{\mathit{\gamma }}_{{\mathrm{\Delta }}_i,m}{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}^{\mathbf{T}}{\mathbf{\Omega }}_{k,m}{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}^{\mathbf{T}}\\ {\displaystyle }& {\displaystyle }={\sum }_{m=\mathrm{1},\mathrm{\dots },M_a}{\mathit{\gamma }}_{\mathrm{a},m}(\mathrm{\Delta }{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}^{\ast }{)}^{\mathbf{T}}{\mathbf{W}}_{i,m}^{-\mathrm{1}}\mathrm{\Delta }{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}^{\ast }\\ {\displaystyle }& {\displaystyle }={\sum }_{n=\mathrm{1},\mathrm{\dots },N_f}{\mathit{\gamma }}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}{\mathrm{\Psi }}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{m=\mathrm{1},\mathrm{\dots },\mathrm{4}}{\mathit{\gamma }}_{{\mathrm{\Delta }}_i,m}{\mathrm{\Psi }}_{{\mathrm{\Delta }}_{\mathit{i}\mathbf{,}\mathit{m}}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{m=\mathrm{1},\mathrm{\dots },\mathrm{6}}{\mathit{\gamma }}_{a_i,m}{\mathrm{\Psi }}_{{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}}\end{array}\end{array}$$

Thus, based on Eqs. (A4) and (A12), the multi-pixel solution for developed multi-platform retrieval is the value minimizing the following quadratic form:

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\displaystyle }\mathrm{\Psi }\left(\mathit{a}\right)& {\displaystyle }={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}{\mathrm{\Psi }}_i\left({\mathit{a}}_i\right)+{\mathrm{\Psi }}_{\text{inter-pixel}}\\ {\displaystyle }& {\displaystyle }={\sum }_{i=\mathrm{1},\mathrm{\dots },N_{\text{pixel}}}[{\sum }_{n=\mathrm{1},\mathrm{\dots },N_f}{\mathit{\gamma }}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}{\mathrm{\Psi }}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{m=\mathrm{1},\mathrm{\dots },\mathrm{4}}{\mathit{\gamma }}_{{\mathrm{\Delta }}_i,m}{\mathrm{\Psi }}_{{\mathrm{\Delta }}_{\mathit{i}\mathbf{,}\mathit{m}}}\\ {\displaystyle }& {\displaystyle }+{\sum }_{m=\mathrm{1},\mathrm{\dots },\mathrm{6}}{\mathit{\gamma }}_{a_i,m}{\mathrm{\Psi }}_{{\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}}]\\ {\displaystyle }& {\displaystyle }+{\sum }_{k=t,x,y}{\mathit{a}}^{\mathbf{T}}{\mathit{\gamma }}_k{\mathbf{\Omega }}_k\mathit{a}\sim min\end{array}\end{array}$$

As a result, the developed retrieval can be controlled by the following Lagrange multipliers (“weights”):

• N_f weights ${\mathit{\gamma }}_{{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{n}}}$ determining the contribution of different instruments or different channels (those parameters are considered pixel-independent);

• up to 4 weights ${\mathit{\gamma }}_{{\mathrm{\Delta }}_i,m}$ determining the contribution of in-pixel smoothness constraints for the vectors ${\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}$ representing some continuous spectral functions or mixtures a_comp

• up to 6 weights ${\mathit{\gamma }}_{a_i,m}$ determining the contribution of in-pixel a priori estimates ${\mathit{a}}_{\mathit{i}\mathbf{,}\mathit{m}}^{\ast }$ for vectors ${\mathit{a}}_{C_{\mathrm{v}}},{\mathit{a}}_{\text{comp}},a_{\mathrm{h}},{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}$. It should be noted here that usually only one a priori in-pixel limitation (either smoothness constraints or a priori estimates) is used for the same parameter.

• up to 6 weights ${\mathit{\gamma }}_{x,m}={\mathit{\gamma }}_{y,m}$ determining the contribution of a priori inter-pixel smoothness constraints on x,y horizontal variability of vectors ${\mathit{a}}_{C_{\mathrm{v}}},{\mathit{a}}_{\text{comp}},a_{\mathrm{h}},{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}$;

• up to 6 weights γ_t,m determining the contribution of a priori inter-pixel smoothness constraints on t-temporal variability of vectors ${\mathit{a}}_{C_{\mathrm{v}}},{\mathit{a}}_{\text{comp}},a_{\mathrm{h}},{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}$;

It should be noted that all Lagrange parameters γ_i can be quantitatively determined as ${\mathit{\gamma }}_i={\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}/{\mathit{\sigma }}_i^{\mathrm{2}}$ (Eq. A6) from knowledge of measurements or a priori data covariance matrices. However, in practice, the knowledge of those covariance matrices is usually quite uncertain. Additionally, some assumptions of the methodology, e.g., dominance of the measurement errors over uncertainties of the used forward model, etc. (see discussion by Dubovik et al., 2021a). Therefore, in practical application of the methodology presented here, the “weights” γ_i are usually adjusted using sensitivity studies and tuning based on acquired experience in retrieval tests.