## 3 Materials and Methods
### Emulation  
Instead of entering the computationally expensive coupled PROSAIL-MODTRAN into GSA, we used an emulated version of these coupled models. Emulation is a statistical learning technique used to estimate model simulations when the model under investigation is too computationally costly to be run many times [71]. The basic idea is that an emulator uses a limited number of simulator runs, i.e., input-output pairs (corresponding to training samples), to train a machine learning regression algorithm (MLRA) in order to infer the values of the complex simulator output given a yet-unseen input configuration. These training data pairs should ideally cover the multidimensional input parameter space using a space-filling sampling algorithm, e.g., LHS. Once the emulator is built, it is not necessary to perform any additional runs of the model; the emulator computes the output that is otherwise generated by the RTM. When it comes to emulating RTM spectral outputs, however, the challenge lies in delivering a full spectrum. This implies that the MLRA should be able to generate multiple outputs to reconstruct a full spectral profile, which is not a trivial task. For instance, the contiguous spectral profile between 400 and 2500 nm consists of over 2000 bands when binned to 1 nm resolution. Only some MLRAs can obtain multi-output models, but that typically lead to highly complex models with long training time and certain risk of overfitting because of model over-representation, e.g., as with neural networks. A workaround solution was developed that enables the regression algorithms to cope with large spectroscopy datasets by taking advantage of the so-called _curse of spectral redundancy_, i.e., the Hughes phenomenon [72]. Since spectroscopy data usually shows a great deal of collinearity, it implies that such data can be compressed to a lower-dimensional space through dimensionality reduction techniques such as principal component analysis (PCA) [73]. Accordingly, spectroscopy data can be converted into components, which are only a fraction of the original amount of bands, and implies that the multi-output problem is greatly reduced to a number of components that preserve the spectral information content (see also [74; 75; 76; 77]). Afterwards, the components are then reconstructed again to spectral data [31; 32; 33; 36; 77]. In earlier RTM emulation evaluation studies [31; 32; 33], various MLRAs were analyzed on their predictive performance. In each of these studies Gaussian processes regression (GPR) [78] was evaluated as the top performing one. Although its superior performance went somewhat at the expense of processing speed as opposed to other MLRAs, it runs numerous times faster than the original RTM [31; 32; 33]. GPR is a probabilistic kernel method, and has been widely used for retrieval of biogeophysical variables and emulation applications [79; 80; 81]. Kernel methods in machine learning owe their name to the use of kernel functions [82; 83; 84]. These functions quantify similarities between input samples of a dataset. Similarity reproduces a linear dot product (scalar) computed in a possibly higher dimensional feature space, yet without ever computing the data location in the feature space. GPR generalize Gaussian probability distributions in function spaces [78]. The prediction and the predictive variance of the model for new samples are given by:  
\[\widehat{\mathbf{f}}(\mathbf{x}_{q}) =\sum_{i=1}^{n}{k(\mathbf{x}_{i},\mathbf{x}_{q})} \tag{10}\] \[\mathbb{V}[\widehat{\mathbf{f}}(\mathbf{x}_{q})] =k(\mathbf{x}_{q},\mathbf{x}_{q})-_{*}^{\top}{(K+\sigma_{n}^{2} \mathbf{I})_{*}^{-1}} \tag{11}\]  
where \(k(\cdot,\cdot)\) is a covariance (or kernel function), \({}_{*}\) is the vector of covariances between the query point, \(\mathbf{x}_{q}\), and the \(n\) or training points, and \(\sigma_{n}^{2}\) accounts for the noise in the training samples.