Sentinel-2 Surface Reflectance Imagery

• Sentinel-2 surface reflectance imagery is atmospherically corrected, providing multi-spectral measurements from 10 m to 60 m for detailed terrestrial and aquatic analysis.
• It employs advanced processing methods like orthorectification, cloud masking, BRDF-harmonization, and deep learning-based gap-filling to produce consistent bottom-of-atmosphere reflectance products.
• This imagery supports applications such as crop monitoring, hydrologic modeling, and ecosystem assessment, with robust metrics (e.g., R² up to 0.89) validating its effectiveness.

Sentinel-2 surface reflectance (SR) imagery comprises multi-spectral, atmospherically corrected observations of the Earth’s surface, acquired by the European Space Agency's Copernicus Sentinel-2 A/B satellites. These Level-2A products represent bottom-of-atmosphere reflectance at spatial resolutions ranging from 10 m to 60 m, enabling fine-scale, multi-temporal analysis of terrestrial and aquatic environments. Underlying processing chains routinely include orthorectification, atmospheric correction, cloud masking, geometric harmonization, and, increasingly, temporal gap-filling and super-resolution. Such standardized SR datasets serve as the substrate for advanced modeling of physical, biological, and ecological dynamics through both statistical and deep learning approaches.

1. Sentinel-2 Surface Reflectance Data Structures

Sentinel-2 Level-2A products encode SR for ten principal spectral bands (excluding aerosol and cirrus channels), with native resolutions summarized as:

Band	ID	Resolution	λ (nm)
Blue	B02	10 m	496.6
Green	B03	10 m	560.0
Red	B04	10 m	664.5
Red Edge 1	B05	20 m	703.9
Red Edge 2	B06	20 m	740.2
Red Edge 3	B07	20 m	782.5
NIR	B08	10 m	835.1
SWIR 1	B11	20 m	1613.7
SWIR 2	B12	20 m	2202.4

All bands are typically resampled to a common 10 m grid via bicubic interpolation when full multispectral stacks or time series are required for analysis (Frion et al., 2023). Level-2A BOA SR is computed from Top-of-Atmosphere (TOA) radiances using processors such as ESA Sen2Cor, which correct for atmospheric path radiance, aerosol loading, water vapor, ozone, terrain effects, and adjacency effects (Revillion et al., 2024). Metadata includes acquisition geometry (sun/view angles), processing baseline, scene classification, and quality masks, facilitating subsequent harmonization and masking operations.

2. Geometric and Atmospheric Correction Workflows

The conversion from TOA to BOA reflectance is a non-trivial radiative transfer problem. Sen2Cor implements:

$$ρ_{\text{surf}}(λ) = \frac{ρ_{\text{toa}}(λ) - ρ_{\text{path}}(λ)}{T_{\text{down}}(λ)×T_{\text{up}}(λ)}$$

where $ρ_{\text{path}}$ is path radiance, and $T_{\text{down}}$, $T_{\text{up}}$ are transmittances (Revillion et al., 2024, Martínez-Ibarra et al., 10 Oct 2025). Scene classification maps, cloud probability masks, and additional layers (water vapor, aerosol optical thickness) are generated. Reprojection, if necessary, uses GDAL/PROJ Transverse Mercator routines for precise UTM/WGS84 alignment.

Atmospheric correction in aquatic systems frequently employs neural network-based modules, such as C2RCC, to retrieve water-leaving reflectance ($ρ_w(λ)$), utilizing both uncorrected TOA and angular-dependent corrections. C2RCC outputs rhow (angular-dependent water leaving reflectance) and rhown (normalized to visible bands) (Martínez-Ibarra et al., 10 Oct 2025).

3. BRDF-Harmonization and NBAR Computation

Surface reflectance is subject to bidirectional anisotropy, necessitating normalization to a reference geometry—Nadir BRDF-Adjusted Reflectance (NBAR). The $c$-factor method achieves this harmonization via a semi-empirical MODIS-derived BRDF model:

$$f_r(\theta_i, \theta_r, \phi) = f_{\text{iso}} + f_{\text{vol}} K_{\text{vol}}(\theta_i, \theta_r, \phi) + f_{\text{geo}} K_{\text{geo}}(\theta_i, \theta_r, \phi)$$

Where $K_{\text{vol}}$ and $K_{\text{geo}}$ are Ross–Li volumetric and geometric kernels. The $c$-factor computes:

$$c(\theta_i, \theta_r, \phi) = \frac{f_r(\theta_i,0,\phi)}{f_r(\theta_i,\theta_r,\phi)}$$

And therefore,

$$\text{NBAR} = c(\theta_i, \theta_r, \phi) \cdot ρ(\theta_i, \theta_r, \phi)$$

Open-source pipelines such as sen2nbar provide high-level functions for batch harmonization of SAFE files and Earth System Data Cubes (ESDCs), supporting both per-band GeoTIFF/COG/NetCDF outputs and lazy, parallelized xarray workflow integration (Montero et al., 2024).

4. Cloud Masking, Index Extraction, and Automated Chain Integration

Successful analysis of SR imagery requires robust cloud identification. Scene classification (SCL) and cloud probability masks (MSK_CLDPRB) from Sen2Cor are thresholded and morphologically processed to yield binary cloud masks, with multiple variants tailored for land-cover-specific reliability (Revillion et al., 2024). Non-cloud pixels serve as input for spectral index computation, with well-established indices (NDVI, NDWI, EVI, NBR, MNDWI) derived via canonical formulas:

• $$NDVI = (\rho_{NIR} - \rho_{R}) / (\rho_{NIR} + \rho_{R})$$
• $$EVI = G\,(\rho_{NIR} - \rho_{R}) / (\rho_{NIR} + C_1 \rho_{R} - C_2 \rho_{B} + L)$$

Sen2Chain orchestrates end-to-end, parallelized operations from download, L2A correction, cloud masking, index extraction, up to time-series export, supporting large-area, high-frequency environmental monitoring (Revillion et al., 2024).

5. Temporal Gap-Filling, Forecasting, and Advanced Modeling

Operational SR imagery is subject to temporal gaps, primarily due to cloud cover. Recent work employs deep models to learn reflectance dynamics for gap-filling, forecasting, and assimilation. Koopman operator-based autoencoders use pixel-wise multispectral time series, embed them in a latent space, and learn a linear latent propagation operator $K$:

• For $x_{i,t} ∈ ℝ^{10}$, measurement $$y_{i,t} = [x_{i,t+1}; x_{i,t+1} - x_{i,t}] ∈ ℝ^{20}$$,
• Encoder $\phi$, decoder $\psi$,
• $z_{t+1} = K z_t$, and prediction $y_{i,t+τ} ≈ \psi(K^{τ} \phi(y_{i,t}))$ (Frion et al., 2023).

Curriculum loss functions (combining short- and long-term objectives plus orthogonality regularization on $K$) support stable, long-horizon prediction. Variational assimilation and gradient-based fitting yield sub-$10^{-3}$ MSE on multi-year time series, outperforming classic interpolators (Cressman shows $>25\%$ higher error).

Complementary approaches use temporal sequence models—attention-augmented BiLSTM architectures—to synthesize missing imagery. By leveraging windows of clear-sky acquisitions and per-date time-difference features, these networks minimize RMSE/MAE/MAPE; direct “gap-free” SR imagery generation for any requested date has been operationalized (Zhao et al., 2024).

6. Super-Resolution: Cross-Band Geometry Models

SR bands in Sentinel-2 imagery natively consist of mixed spatial resolutions, complicating fine-scale analysis. The "band-independent geometry" model separates pixel geometry from band-dependent reflectance, fitting shared spatial weights ($g_i$) across high-res bands and then unmixing/sharpening lower-res bands (20 m, 60 m) to produce a full stack of 10 m reflectance layers. A local linear predictor $V$ and geometric mean ratio sharpening combine with band registration and optional rescaling to guarantee conservation of coarse-resolution averages (Brodu, 2016). C++/Python implementations enable scalable processing; performance metrics (Q, ERGAS, SAM) demonstrate high fidelity to original imagery across diverse land covers.

7. Application Domains and Best Practices

Sentinel-2 SR imagery underpins advanced applications including multitemporal crop monitoring, hydrologic and canopy process modeling, and aquatic ecosystem assessment. In predictive mapping of water column chlorophyll-a, C2RCC-processed reflectances, robust cloud masking, stratified spatial aggregations, and ML/DL model ensembles (RF, XGBoost, CatBoost, MLPs, KNN) yield high R² (surface: up to 0.89, deeper layers: up to 0.81) and accurately reconstruct event-scale dynamics (Martínez-Ibarra et al., 10 Oct 2025).

Best practices encompass:

• Use atmospherically corrected L2A products for all critical analyses.
• Harmonize geometry via NBAR before computation of temporal indices or learning models.
• Apply spatial aggregations (e.g., 5×5–15×15 uniform mean) where pixel noise is problematic.
• Systematically compare gap-filling and forecasting models against linear and persistence baselines.
• Version-control pipeline configurations, processing environments, and job definitions for reproducibility.

This synthesis of acquisition protocols, geometric/atmospheric processing, harmonization techniques, time-series modeling, gap-filling strategies, and super-resolution algorithms defines the state-of-the-art in Sentinel-2 surface reflectance imagery, ensuring robust, scalable, and physically consistent derivatives for scientific, operational, and policy use cases.