## 3 Materials and Methods
### Global Sensitivity Analysis (GSA)  
In order to identify the driving vegetation and atmospheric variables having an impact on TOA radiance, we first conducted a global sensitivity analysis (GSA) of the TOA radiance simulations from the _analysis_ dataset. Most GSA methods are variance-based methods, which decomposes the variance of the model output into fractions that can be attributed to inputs or sets of inputs [28; 65]. While the [PERSON]' method [66] pioneered in developing a variance-based GSA method, a modified version was proposed by [67], which proved to be effective in identifying the so-called [PERSON]'s  
Figure 4: Typical error in surface reflectance after Sen2 Cor atmospheric correction algorithm. See Table VIII in [50] for detailed information. sensitivity indices. These indices quantify both the main sensitivity effects (first-order effects: \(S_{i}\), i.e., the contribution to the variance of the model output by each input variables) and total sensitivity effects (\(S_{Ti}\), i.e., the first-order effect plus interactions with other input variables) of input variables. This method has been applied here. A description according to [68] is given below. Formally, we have a model \(y=f(\mathbf{x})\), where \(y\) is the model output, and \(\mathbf{x}=[x_{1},x_{2},...,x_{k}]^{\top}\) is the input feature vector. A variance decomposition of \(f(\cdot)\) as suggested by [PERSON] [66] is:  
\[\mathbb{V}(y)=\sum_{i=1}^{k}V_{i}+\sum_{i=1}^{k}\sum_{j=i+1}^{k}V_{ij}...+V_{1,...,k}, \tag{4}\]  
where \(\mathbf{x}\) is rescaled to a \(k\)-dimensional unit hypercube \(\Omega^{k}\), \(\Omega^{k}\) = \(\{\mathbf{x}|0\leq x_{i}\leq 1,i=1,\ldots,k\}\); \(\mathbb{V}(y)\) is the total unconditional variance; \(V_{i}\) is the partial variance or'main effect' of \(x_{i}\) on \(y\) and given by the variance of the conditional expectation \(V_{i}=\mathbb{V}[\mathbb{E}(y|x_{i})]\); \(V_{ij}\) is the joint impact of \(x_{i}\) and \(x_{j}\) on the total variance minus their first-order effects. Here, the first-order sensitivity index \(S_{i}\) and total effect sensitivity index \(S_{Ti}\) are given as [28]:  
\[S_{i}=\frac{V_{i}}{\mathbb{V}(y)}=\frac{\mathbb{V}[\mathbb{E}(y|x_{i})]}{ \mathbb{V}(y)} \tag{5}\]  
and:  
\[S_{Ti}=S_{i}+\sum_{j\
eq i}S_{ij}+...=\frac{\mathbb{E}[\mathbb{V}(y|x_{\sim i })]}{\mathbb{V}(y)}, \tag{6}\]  
where \(x_{\sim i}\) denotes variation in all input variables and \(x_{i}\), \(S_{ij}\) is the contribution to the total variance by the interactions between variables. Following [PERSON] et al. (2010) [67], to compute \(S_{i}\) and \(S_{Ti}\), two independent input variable sampling matrices \(\mathbf{P}\) and \(\mathbf{Q}\) of dimensions \(N\times k\) are created, where \(N\) is the sample size and \(k\) is the number of input variables. Each row in matrices \(\mathbf{P}\) and \(\mathbf{Q}\) represents a possible value of \(\mathbf{x}\). The variable ranges in the matrices are scaled between 0 and 1. The Monte Carlo approximations for \(\mathbb{V}(y)\), \(S_{i}\) and \(S_{Ti}\) are defined as follows [67; 69]:  
\[\hat{\mathbb{V}}(y)=\frac{1}{N}\sum_{j=1}^{N}(f(\mathbf{P})_{j})^{2}-\hat{f}_ {0}^{2},\quad\hat{f}_{0}=\frac{1}{N}\sum_{j=1}^{N}f(\mathbf{P})_{j}, \tag{7}\]  
and:  
\[\widehat{S}_{i}=\frac{1}{N}\sum_{j=1}^{N}\frac{f(\mathbf{Q})_{j}(f(\mathbf{P} _{\mathbf{Q}}^{(i)})_{j}-f(\mathbf{P})_{j})}{\hat{\mathbb{V}}(y)} \tag{8}\]  
and:  
\[\widehat{S_{Ti}}=\frac{1}{2N}\sum_{j=1}^{N}\frac{(f(\mathbf{P})_{j}-f(\mathbf{ P}_{\mathbf{Q}}^{(i)})_{j})^{2}}{\hat{\mathbb{V}}(y)}, \tag{9}\]  
where \(\widehat{\cdots}\) is the estimate; \(\hat{f}_{0}\) is the estimated value of the model's output; we abused notation by defining \(f(\mathbf{P})\) as all outputs for row vectors in \(\mathbf{P}\); \(\mathbf{P}_{\mathbf{Q}}^{(i)}\) represents all columns from \(\mathbf{P}\) except the \(i^{th}\) column which is from \(\mathbf{Q}\), using a radial sampling scheme [70].