## 3 Materials and Methods
### Hybrid Retrieval Schemes  
When it comes to selecting a biophysical variable retrieval method for processing large images such as Sentinel-2 (S2), it requires models that are fast, robust and easily applicable. Based on a systematic comparison of parametric, non-parametric and RTM-inversion retrieval methods taking both accuracies and run-time into account [86], it was concluded that hybrid retrieval schemes, i.e., machine learning methods trained by RTM simulations, can achieve both accurate and fast estimates. Regarding the used MLRA, similar as in emulation, GPR was evaluated as a powerful method for mapping applications [38; 39; 87]. Starting with Equation (10), we used a scaled Gaussian kernel function:  
\[k(\mathbf{x}_{i},\mathbf{x}_{j})=\
u\exp\bigg{(}-\sum_{b=1}^{B}\frac{\left( \mathbf{x}_{i}^{(b)}-\mathbf{x}_{j}^{(b)}\right)^{2}}{2\sigma_{b}^{2}}\bigg{)} +\delta_{ij}\cdot\sigma_{n}^{2}, \tag{14}\]  
Regarding retrieval, three important properties of the method are worth stressing here. First, the obtained weights \(w\) after optimization gives the relevance of each spectrum \(\mathbf{x}_{i}\) (see [88] for extended equations). The predictive mean is essentially a weighted average of the vegetation biophysical parameter values associated with the training samples closest to the test sample. Second, the inverse of \(\sigma_{b}\) represents the relevance of band \(b\). Intuitively, high values of \(\sigma_{b}\) mean that relations largely extend along that band hence suggesting a lower informative content. These features have been extensively studied in [87; 88] and proved to be valuable for gaining insight into relevant bands. Third, and particularly of interest for mapping applications, a GPR model provides not only a per-pixel prediction, but also an uncertainty (or confidence) level for the prediction. Hence, uncertainty intervals are directly delivered along with the variable estimates, which enables to assess the model transferability in space and time [86; 88].  
We assessed the performance on biophysical variable retrieval on the three reference scenarios previously described in Section 3.2. The MLRA retrieval toolbox was first used to train and to validate an MRLA from the _retrieval_ datasets and then to apply the trained model to retrieve biophysical variables from the _reference_ datasets. Based on experience from earlier retrieval exercises [32; 33], the two _retrieval_ databases were split into 70% for the training and 30% for the validation of the GPR retrieval algorithms. The _retrieval_toc_ dataset was used to train and validate one GPR model for the retrieval of biophysical variables from surface reflectance. The _retrieval_toa_ was instead used for the retrieval from TOA radiance. Goodness-of-fit statistics were calculated to assess the GPR models' capability to retrieve accurately biophysical variables. The error difference between the reference and the retrieved biophysical variables is calculated for each of the 5000 samples in the _reference_ database and the histogram plotted. In addition, Pearson's correlation coefficient (\(R^{2}\)) and root-mean-square error (RMSE) are calculated.