Bayesian inversion of a forest reflectance model using Sentinel-2 and Landsat 8 satellite images

The inversion of reﬂectance models is a generalizable tool to obtain estimates on forest biophysical pa- rameters, such as leaf area index, with theoretically little information need from a study area, instead relying on the knowledge about physical processes in the forest radiation regime. The use of prior infor- mation can greatly improve the reﬂectance model inversion, however, the literature does not yet provide much information on the selection of priors and their inﬂuence on the inversion results. In this study, we used a Bayesian approach to invert the PARAS forest reﬂectance model and retrieve leaf area index from Sentinel-2 MSI and Landsat 8 OLI multispectral satellite images. The PARAS model is based on the the- ory of spectral invariants, which describes the inﬂuence of wavelength-independent parameters on forest radiative transfer. The Bayesian inversion approach is highly ﬂexible, provides uncertainty quantiﬁcation, and enables the explicit incorporation of prior knowledge into the inversion process. We found that the choice of prior information is crucial in inverting a forest reﬂectance model to predict leaf area index. Regularizing and informative priors for leaf area index strongly improved the predictions, relative to an uninformative prior, in that they counteracted the saturation effect of the optical signal occuring at high values for leaf area index. The predictions of leaf area index were more accurate for Landsat 8 than for Sentinel-2, due to potential inconsistencies in the visible bands of Sentinel-2 in our data, and the higher spectral resolution.


Introduction
Forest reflectance models aim at providing a generalized tool to researchers to both understand and simulate the reflectance of a forest stand, given a set of forest biophysical parameters, for example the leaf area index (LAI), that influence the forest radiation regime. These models are typically ill-posed, which makes their inversion, i.e., the retrieval of forest biophysical variables, a non-trivial problem. A number of inversion methods have been suggested in the literature [e.g. 19 ], and most methods produce point estimates as output of the forest reflectance model inversion. These point estimates are obtained by applying one of several radiative transfer inversion methods, that are, look-up tables, optimization methods, and neural networks, among others. These inversion methods have been proven to provide reasonable retrieval accuracy of forest biophysical variables, given that the construction of the look-up table, or the training data for a neural network, cover the ecologically reasonable parameter range of radiative transfer model inputs [19] . However, the accuracy with which many inversion methods predict forest biophysical parameters strongly depends on the way they are parameterized, or what initial guess one chooses. Optimization algorithms in particular may be challenging due to their need for an initial value guess, which can lead to a significant error [45] . Recent advances in computations and algorithms related to Bayesian inference have made this inversion approach more accessible. The Bayesian inversion approach leverages the theoretical foundation of Bayes' theorem and usually Markov Chain Monte Carlo (MCMC) algorithms to estimate the probability distribution of unknown parameters in a reflectance model. In the case of reflectance model inversion, those unknowns are the forest structural and spectral parameters, conditional on observed bottom-of-atmosphere (BoA) reflectance. The Bayesian framework features, among others, two properties that make it particularly suitable for reflectance model inversions, a point estimate, which provides inherent uncertainty quantification, and (2) the possibility to explicitly incorporate existing knowledge on parameters into the model itself. A few studies have applied a Bayesian inversion approach to forest reflectance models [45,46,51] or the leaf optical model PROSPECT [37] with promising results, particularly in terms of uncertainty quantification. Varvia et al. [45] has laid out the ground work for a Bayesian inversion approach using the PARAS forest reflectance model [34] .
The PARAS model is based on the theory of spectral invariants, which describes the influence of canopy structure on the forest radiation regime. The model parameterizes the forest radiation regime using only few spectral and structural parameters, which makes it particularly suitable for a Bayesian inversion. Furthermore, PARAS has already shown promising results in the inversion of hyperspectral data to retrieve leaf area index, and leaf and understory optical properties [45,46] .
Physically-based reflectance model inversion usually sets out to retrieve forest biophysical parameters, such as leaf area index or other canopy structural or biochemical parameters, from observed airborne or spaceborne imagery. This biophysical retrieval approach is only in part competing with empirical, or statistical, retrieval approaches because their relative advantages are rather different. The physical approach is more suitable for large area, multi-sensor and low field data scenarios, whereas the empirical approach is more suitable for regional applications with high prediction accuracy requirements [38] . A key property of physical models is their attempt to be sensor independent, which is becoming more important with every new satellite mission to be able to maintain a consistent theoretical foundation for combining data from different sensors.
In optical remote sensing of vegetation, two recent satellite missions have attracted a lot of attention in science and industry communities; the high resolution optical satellite sensors onboard Sentinel-2a and -2b (Multi Spectral Instrument, MSI), and onboard Landsat 8 (Operational Land Imager, OLI). Hereafter, we use the names of the satellite platforms and their respective sensors synonymously, abbreviated as S2 and L8. These two multispectral sensors share bands of similar spectral sensitivity, which appear to generally agree well [e.g. [4,22,25] ], despite some studies reporting differences particularly in the red band, where S2 appears to exhibit lower reflectance in some cases [20,44,47] . S2 further features three bands in the red edge spectral region, which has been shown to slightly improve retrievals of forest biophysical parameters [20] . Besides the novelty of the additional red edge bands of S2, both sensors have an unprecedented spatial resolution for freely available satellite imagery; 10/20m resolution for vegetation bands of S2, and 30m resolution for L8. From a Bayesian reflectance model inversion perspective, these two sensors are interesting because previous studies used PARAS model only on hyperspectral data [46] . The effect of spectral resolution was previously studied by Shiklomanov et al. [37] at leaf level, where they pointed out that inversion uncertainty increases with decreasing spectral resolution. They also noted that L8 in particular has, despite only subtle differences in spectral response, performed better in the model inversion than its predecessors Landsat 5 and Landsat 7, potentially hinting that the multispectral resolution of L8 to be in a particularly interesting range.
Despite many effort s f or sensor-intercomparison in f orest applications in general [e.g. 1,12 ], studies that analyze the relative potential of S2 and L8 for forest biophysical parameter retrieval remain few, and are limited to empirical modeling approaches [e.g. 20 ]. To our knowledge, this is the first study to apply a physicallybased forest biophysical parameter retrieval using S2 and L8 data. Besides the attempt to compare S2 and L8 in their potential for forest reflectance model inversion, we further set out (1) to determine whether this inversion approach is feasible at this level of spectral resolution, (2) to estimate the inversion quality through the accuracy of leaf area index retrievals, and (3) to analyze the influence of the choice of prior distribution on the inversion accuracy.

Study area, field and lidar plots
The study area in Suonenjoki, Central Finland (45 • 18'N, 123 • 21'W) is dominated by the typical Finnish boreal forest tree species Norway spruce ( Picea abies (L.) Karst.), Scots pine ( Pinus sylvestris L.), and Silver and Downy birch ( Betula pendula Roth and Betula pubescens Ehrh.). The area is characterized by diverse forest structure and density, with average stand size of about 2 ha. Our reference data consisted of 746 lidar plots for which the leaf area index ( L ) was predicted using airborne lidar data. The LiDAR data were collected on 4 September, 2014 using a Leica ALS70-HA laser scanning system. It has a wavelength of 1064 nm and a beam divergence of 0.22 mrad. The flight altitude was 20 0 0 m above ground level, and the scan angle was constrained to ± 15 degrees. The nominal pulse density of the lidar data was 0.75 pulses per square meter, with a footprint size of about 44 cm. There were 19 plots with field-measured L within the scanned area which were used to train regression models for predicting L for the lidar plots. L was retrieved from the fraction of all echos over 1.3 m above ground, and canopy cover from the fraction of all first and single echos over 1.3 m above ground, based on Korhonen et al. [21] . The lidar plots had a square shape with edge length of 60 m, and were placed over a set of nine (S2) or four (L8) pixels. Due to geolocation differences of the images, the plot center points were located around 16 m away from each other. More details are available in Korhonen et al. [20] , where the dataset was first described and is available from their supplemental materials.
The S2 image was obtained from Sentinel Scientific Data Hub [10] as Level-1C product (geometrically corrected top-ofatmosphere reflectance). An atmospheric correction was done using the Sen2cor module version 2.2.1 within the Sentinel 2 Toolbox (S2TBX), Sentinel Application Platform (SNAP) version 4.0.2. In this study, spectral bands with 10 m and 20 m resolution were used. S2's band 8 was omitted from the analysis due to potential inconsistencies with L8, as recommended by Li et al. [22] . Geometrically corrected Landsat 8 Surface (bottom-of-atmosphere, BoA) Reflectance Climate Data Records [27] were downloaded from the United States Geological Survey [43] . The square-shaped lidar plots contained the average of nine (S2) and four (L8) pixels.

PARAS model and parameterization
The forest reflectance model PARAS [34] is based on the concept of spectral invariants, specifically the photon recollision probability, and relies only on a handful of forest structural and spectral parameters. The model simplicity makes PARAS suitable for inversion with Markov Chain Monte Carlo, as the ill-posedness of a radiative transfer model is expected to increase with the number of parameters, as parameter identification becomes more problematic [e.g. 32 ]. PARAS models forest reflectance as or in other words, the forest stand bidirectional reflectance factor ( r ) consists of a reflectance contribution from the ground or understory ( r G ) and the tree canopy ( r C ). Both contributions depend on the illumination and viewing zenith angles ( θ i and θ v , respectively), and the wavelength ( λ).
The understory contribution describes the proportion of sunlit ground area that is visible to the sensor: with T , canopy transmission (also called canopy gap fraction), forming the bidirectional gap probability T ( θ i ) T ( θ v ) in illumination and viewing direction, and the understory hemisphericaldirectional reflectance factor ( ρ G ).
The tree canopy contribution to the reflected signal is modeled as the fraction of intercepted sunlight that is scattered into the viewing direction: with the canopy interceptance i 0 (θ i ) = 1 − T (θ i ) , the canopy albedo ω C and the ratio of backscattered radiation Q . The canopy albedo is based on the photon recollision probability p , or the probability that a photon, upon an interaction with a canopy element, will interact with another element within the canopy. Expressing this assumption mathematically, and using ω L to denote the canopy element albedo, one gets the geo- where k denotes the scattering order (i.e. k = 1 is first-order scattering, or photons that escape the canopy after one interaction). This geometric series converges to with the common ratio ω L p , as the photon has additional interactions with canopy elements and is scattered, and the first term of the series is ω L (1 − p) , which, from a causal perspective, are actually the first and last parts of the sequence. That is, the photon is scattered from its initial interaction within the canopy ( ω L ) and is ultimately escaping the canopy ( 1 − p). Q is the directional parameter describing the fraction of canopyscattered radiation that escapes in an upward direction, i.e. towards the sensor. This backscatter fraction was parameterised as proposed by Mõttus and Stenberg [28] : This parameterisation of Q incorporates an additional spectral invariant parameter, q , which was empirically determined by Mõttus and Stenberg [28] as q = exp (−0 . 1684 L e ) . The effective leaf area index L e in this context is related to the (true) leaf area index as L e = βL via the clumping index β. We use the term effective leaf area index as defined in Black et al. [3] . The photon recollision probability p contains the most information on forest structure in PARAS model, and it depends on both clumping and the angular dependency of canopy gaps over the entire hemisphere. We parameterized p following Stenberg [40] : with the diffuse interceptance i D , the clumping index β and the effective leaf area index L e . The diffuse interceptance in PARAS is described as The angular relationship between L e and T was modeled using the modified Beer's law with a constant clumping index: with the projection of unit leaf area into a zenith angle G ( θ ). Parameterizing G requires estimates of the leaf angle distribution, which we incorporated in our inversion model. G depends on leaf angle ( θ L ) as [49] A in this equation describes the projection area ratio of leaves with zenith angle θ L into the direction θ , and incorporates assumptions on the shape of leaves, and f is the frequency distribution of leaf angles. We modeled A separately for deciduous and conifer trees. For deciduous leaves, we used the projection function for flat leaves [48] For conifer leaves, we used the function for non-flat leaves proposed by Stenberg [39] : (11) We modeled f also separately for deciduous and conifer leaves. Deciduous trees in Finland are predominantly birch species, for which data was available from Pisek et al. [31] , who also clearly recommended the use of planophile leaf angle distributions for boreal broadleaved tree species. Therefore, we used a two-parameter Beta distribution [13] to model the leaf angle distribution of deciduous leaves from the existing data as: with the Beta function B , and two parameters μ and ν, which depend on the mean and standard deviation of θ L . This distribution was recommended by Wang et al. [48] for cases where data on mean and standard deviation is available. There was a lack of data for boreal conifer tree species leaf angle distributions, hence we assumed a spherical leaf angle distribution, following De Wit [6] f (θ L ) = sin (θ L ) .
We approximated the integrals in Eqs. (9) -(11) by quadrature using 10 evenly spaced locations, and Eqs. (7) and (8) at five discrete angles with associated weights, corresponding to the method used for the LAI-20 0 0 instrument [23] . In a nutshell, the PARAS model has the unknown parameters , q, p . The latter two parameters, q and p , require further parameterization, and in particular p requires a complete parameterization of forest structure, namely through the effective leaf area index, the clumping index, and the leaf angle distribution. We assumed a planophile leaf angle distribution for deciduous trees [31] , and a spherical distribution for conifer trees. The PARAS model does not inherently model tree species, however, it is possible to model tree species through a mixture parameter for leaf optical properties and clumping index. We distinguished only deciduous and conifer tree species to avoid parameter identification problems between the two conifer species, Norway spruce and Scots pine, and refer to the species parameter as conifer share, which is defined as the share of conifer species in the leaf area index. Parameters that are modeled by species are denoted by a subscript "c" for conifer, and "d" for deciduous.

Bayesian reflectance model inversion
The forward problem of PARAS, which simulates forest stand reflectance from given forest structural and spectral parameters, can be stated as r = h (x ) + . In this context, r ∈ R m ×n is the observation matrix, it contains the satellite pixel measurements for the number of plots m = 746 and the number of sensor bands n , i.e. n = 6 for L8, and n = 9 for S2. structural properties L e , β c , β d , p c were modeled for each plot i in the study area. We chose this hierarchical modeling approach to leverage assumed similarities between plots.

The vector of unknown parameters
The following elaboration on the Bayesian inversion of PARAS closely follows Varvia et al. [45] . In the Bayesian approach, both observed and unknown parameters are modeled as random variables. We therefore denote the prior probability density by P ( x ), which contains available information on x before any measurements were done. We relied on readily available information from peer-reviewed literature in combination with some practical constraints to select suitable prior distributions for each parameter in x . More details on prior selection are given in Section 2.5 . In Bayesian inference, the prior density is updated by the information gained from satellite measurements using Bayes theorem where P ( x | r ) is the posterior probability density, P ( r | x ) is the likelihood function, and P ( r ) is the posterior marginal density. The latter term serves to normalize the posterior, and can be neglected in the further analysis. The posterior density is the full solution to the Bayesian inverse problems, i.e. it is a joint probability density of all forest structural and spectral parameters in PARAS. The likelihood function contains the information from the measurement, i.e. the satellite bands' observed BoA reflectance factor for the lidar plots. The likelihood is coincident with the forward problem, i.e. r = h (x ) + , where ∈ R n is an additive error term, containing both model error and measurement noise. In the case of additive error, the likelihood function becomes with the error density P modeled as a multivariate normal distribution with zero mean, and an uncorrelated covariance with standard deviation of 20% of the observations in each band.

Prior choices
To study the effect of prior choices on the inversion result, we chose three different prior distributions for the effective leaf area index: where N denotes a normal distribution, with their respective means and standard deviations. These three priors resemble different modeling approaches. The first approach is uninformative ( Eq. (16) ), it makes no assumption about L e other than imposing a range of possible values. The second prior is informative ( Eq. (17) ), it assumes that the probability distribution for L e in any plot follows the distribution of L e in the geographical region. We approximated this assumption by a normal distribution, based the results of Härkönen et al. [14] , who estimated L e from the Finnish National Forest Inventory data 2004-2008. The third prior is regularizing ( Eq. (17) ), it formulates the probability distribution of L e particularly to regularize the inversion in a way that counteracts a saturation effect of the optical signal at high values of L e . This saturation effect causes the relationship between L e and forest stand reflectance to flatten out in dense forest stands, leading to high uncertainty and potential overestimation of L e in the inversion ( Fig. 3 ). Hence, the third prior was chosen in anticipation and to counteract such a saturation effect.
The prior for conifer clumping index ( β) is based on the assumption that shoot clumping is the main cause for non-random distribution of canopy elements. Following Stenberg et al. [41] , a shoot clumping index can be expressed as four times the silhouette to total area ratio (STAR) of a shoot. The observed range for shoot STAR ranges from about 0.1 to 0.22 [41] , resulting in a clumping index of 0.4 to 0.88. For the prior of β, we chose a mean of 0.6 and a standard deviation of 0.2 to assign the highest probability to the observed range of shoot STAR, but simultaneously allow the clumping index to take values beyond this range, which reflects on the little known influence of higher level clumping. In deciduous trees, shoot clumping is non-existent, and the common assumption is that the random leaf distribution holds. To be able to quantify a potential influence on higher level clumping, e.g. at crown or stand level, we modeled a separate clumping index for deciduous trees centered at 1, with a standard deviation of 0.2. The clumping index for both species were constrained within the range 0.05 and 1.1. The dominant tree species was assumed to be coniferous, expressed through the prior for the conifer share p c with mean of 0.8 and standard deviation of 0.5.
Priors for leaf and understory optical properties were obtained from published measurement data. We used the spectrometer measurements from Hovi et al. [17] for pine, spruce, and birch leaf optical properties. This database contains seasonal measurements of 25 boreal tree species, with reflectance and transmittance measurements reported for several samples and stratified by shoot age, leaf side, and crown position. We selected the measurements for birch on 25 August 2016, and for pine and spruce on 12 September 2016, which were the closest dates to the day of year of the satellite image acquisitions. We added reflectance and transmittance measurements to get leaf albedo, and averaged over shoot age, leaf side, and crown position as well as over individual samples. Finally, the spectra for pine and spruce were averaged to form a single conifer leaf albedo. Understory hemispherical-directional reflectance factor (HDRF) data came from Hovi et al. [16] , where we used the measurements from June 2010, and averaged over the sub-xeric, mesic, and herb rich site types. The spectral data was available in the wavelength region between 400 and 2400 nm. We converted the spectra to the respective sensor bands of S2 and L8 using the sensor-specific spectral response functions. These were obtained from the Sentinel 2 Document Library [9] and Barsi et al. [2] . The resulting mean spectra for leaf and understory were used as the prior means for the deciduous, conifer, and understory optical properties.
Spectral properties were modeled as multivariate normal distributions, centered at their measured values. The covariance between the individual bands was crucial for forest reflectance model inversion to avoid parameter identification problems. In creating a covariance matrix, we followed the correlation concept laid out by Varvia et al. [45] . The correlation matrix R between bands comes from three weighted components: where the unitary matrix means the correlation between all bands, R part means the correlation between selected groups of bands that are assumed similar, and the identity matrix means a random component for each individual band. The above matrices have dimensions n bands × n bands , i.e. 6 × 6 for L8 and 9 × 9 for S2. The factors κ are the weighting coefficients, and we used κ all = 0 . 1 , κ part = 0 . 2 , κ ind. = 0 . 7 . To obtain R part , the sensor bands were grouped, with the visible bands, and the red edge 1 band of S2 forming the first group, the near infrared (NIR) band being its own group, together with red edge 2 and 3 bands of S2, and the shortwave infrared (SWIR) bands 1 and 2 forming the third group. Therefore, for S2, and, for L8, analogously the three partial matrices have dimension 3 × 3, 1, and 2 × 2. The covariance matrix is obtained by = SRS T , where R is the correlation matrix ( Eq. (19) ), and S is the diagonal matrix of standard deviations of individual bands, which we parameterized with 10% of the measured spectral data of each band.

Computations and accuracy assessment
For each of the three models, i.e. the different priors for L e , and per sensor, we sampled 50 0 0 samples in each of 16 parallel chains, with a burn-in phase of 50 0 0 samples. The computations were carried out in Python using the PyMC3 package [36] .
We used a set of different accuracy metrics to assess the accuracy of parameter retrievals from the inversion. The full solution to the inverse problem is the posterior distribution, approximated by MCMC sampling. To describe the posterior, we employed two types of metrics, based on credible intervals or point estimates. We used the credible interval of the posterior marginal distribution that contains 95% of the probability mass for a given parameter. The posterior marginal distribution is approximated by the set of samples only for a single parameter, ignoring the joint probability distribution with other parameters. There is a potentially infinite number of intervals that contain 95% of the probability mass, hence we use the common approach to select the narrowest interval, called the highest posterior density (HPD) interval. The HPD has a straight-forward interpretation, it describes the interval that contains the true value with a probability of 95%. The accuracy metric associated with the HPD is the percentage of plots whose observed values lay within the HPD interval.
To allow comparison with conventional, non-Bayesian accuracy metrics based on point estimates, we used the posterior mode, i.e. the most likely value of the posterior. The conventional accuracy metrics we report are the root-mean square error (RMSE), the mean bias, and the bias corrected RMSE (cRMSE), calculated as follows: where y and ˆ y are the observed and predicted variable, and n is the number of plots. The Python code used for this study is available from the supplementary file (S1).

Results and discussion
We present our findings in two parts. The first is a descriptive analysis of the available data, which is essential in the interpretation of the latter part, which presents the results of our inversion and a discussion on the influence of prior information on the inversion results. Following the presentation and discussion of concrete results is a general discussion on the potential and limitations of our inversion approach.

Descriptive satellite and plot data analysis
The satellite observations for the respective lidar plots of S2 and L8 were mostly in agreement ( Fig. 1 ). However, the BoA reflectance values in the visible bands differ strongly between the sensors. The blue and green bands of S2 show their strongest differences in the darker plots, with the difference decreasing with increasing reflectance. In the red band, the reflectance values of S2 are overall shifted by about 0.01 units, the regression line indicates closer agreement at a reflectance of about 0.1, which is well beyond the highest reflectance observation in the study area. Korhonen et al. [20] already pointed out these reflectance differences between sensors, and noted that the differences likely do not stem from the atmospheric correction procedure. Furthermore, Korhonen et al. [20] pointed out that two other studies [44,47] found similar reflectance differences and speculating at sensor differences as the main reason. Later studies [e.g. 22,50 ] found the differences between the sensors to be much smaller. The ESA Sentinel-2 Team has identified a few product anomalies during the commissioning period from June to December 2015, and the data products were reprocessed twice. The first reprocessing was done using processing baseline 02.01, and, after some product anomalies were found, a second reprocessing was done using baseline 02.04, as documented in the Sentinel-2 Mission Status Report 44, July 2016, available from European Space Agency [11] . We compared the S2 L1C products produced by both processing baselines, and found no difference in the images. Since the sensor differences in visible bands were only found in studies using early Sentinel-2a images (until May 2016), but not in later studies, where the agreement between S2 and L8 seems much better, we suspect that pre-processing of S2 images may have been different during its commissioning phase, or other technical reasons may have been involved. The Bayesian reflectance model inversion approach was able to compensate for the sensor differences in the visible bands to some extent. In the posterior distribution, the leaf albedo was significantly lower for coniferous trees in the visible bands of S2 ( Fig. 2 ). Coniferous stands typically constitute the darkest pixels due to increased clumping [34] , lower reflectance observations may lead to an overestimated conifer share in the inversion, and to a lower leaf albedo of conifers as the dominant species.
The relationship between the leaf area index, and reflectance shows a typical decreasing reflectance with increasing LAI, with a saturation effect at LAI greater than about 3 to 4 ( Fig. 3 ). Similar to the above discussion, the reflectance differs between the sensors in the visible bands. The forest reflectance of stands with small LAI consists mostly of the reflected signal from the bright understory [8,33] . Hence, in accordance with the above discussion on sensor differences, there is little difference in the blue and green band for sparse stands, i.e. plots with LAI close to zero. As the LAI increases, the reflectance of the visible bands decreases and levels out at higher levels for L8 than S2, a difference which is not present in NIR and SWIR bands.

Reflectance model inversion results
The use of prior information in reflectance model inversion in general has been widely accepted [e.g. 19 ]. However, the use of  prior information is usually limited to defining the range of reflectance model parameters that is ecologically meaningful, or using random draws from a field dataset. In terms of Bayesian inference, setting a range is equivalent to assigning a uniform distribution to reflectance model parameters. We suspect that this may not be the best choice, particularly for estimating parameters which exhibit a saturation effect with the remotely sensed signal, such as LAI, because such nonlinearities in the model affect the proper choice of an objective, uninformative prior [e.g. 18 ]. We hypoth-esized that using non-uniform priors improves predictions of optically saturating parameters, that are, parameters which at some range cannot be reliably predicted from the optical signal, or which can only be predicted with immense uncertainty. To test this hypothesis, we employed three reflectance model inversion strategies that rely on different prior parameterizations for LAI. The three strategies, or prior models, were (i) a uniform prior as per common assumption, (ii) a realistic (informative) prior which approximated the LAI of Finnish forests [based on National Forest Inventory Data, 14 ], and (iii) a regularizing prior which was aimed at counteracting the saturation effect of the optical signal at high LAI by favoring small values for LAI.
Overall, both the informative and regularizing prior models performed satisfactory, while the uniform prior was heavily biased ( Fig. 4 ). For L8, the regularizing prior resulted in the best LAI retrievals (bias = 0.03, RMSE = 0.65), compared to the informative prior (bias = −0.19, RMSE = 0.67). While both priors appear very similar in their prediction accuracy of stands with LAI larger than about 2, there is a notable difference in sparse stands with small LAI between 0 and 2. In these sparse stands, the informative prior introduced a visible bias, causing overestimated LAIs especially in the range of LAI values between 0 and 2. The regularizing prior as-signs the highest probability to low LAI values, hence the prediction accuracy is higher for low LAI values, while not compromising predictions of stands with large LAI.
For S2, the informative prior performed slightly better (bias = −0.76, RMSE = 1.13) than the regularizing prior (bias = −0.94, RMSE = 1.34). Here, there is no visible difference for the two priors for LAI values between 0 and 2, as was observed in L8 retrievals. This may be due to the lower reflectance values of S2 in the visible bands, which in turn cause overestimation of LAI in the inversion and have more influence on the results than the difference in LAI priors. The confidence interval metric, the ratio of plots with observed LAI within the HPD, exhibits the same ranking of priors as the point estimate-based metrics in both sensors. In S2, the informative and regularizing priors respectively captured 66% and 63% of the plots correctly. 87% and 85% were the corresponding figures for the regularizing and informative priors for L8, respectively. When considering uncertainty, as the share of plots within the HPDs showed, the difference between the regularizing and informative priors appears less pronounced than with point estimate metrics.
For both sensors, the uninformative priors resulted in significantly worse LAI retrievals (L8 bias = −2.67, RMSE = 2.98; S2 bias = −3.89, RMSE = 4.39). The uninformative case showed a strong influence of the optical signal saturation effect, as the bias increases linearly with increasing LAI (see Fig. 4 ). Both the regularizing and informative priors outperformed the uninformative prior by one or two orders of magnitude in bias. The confidence intervals for LAI in the uniform prior model are relatively wide, however, only a small proportion of plots' LAI values was within the HPDs, i.e. 43% for L8, and only 11% for S2.
The LAI retrieval accuracy for the regularizing and informative priors are satisfactory. The dataset we used in this study was originally used by Korhonen et al. [20] to empirically predict LAI, among other variables, using a plethora of vegetation indices and single and multiple band predictors. They found single bands for both S2 and L8 to produce the best predictions, with the red edge 1 band of S2 predicting LAI with RMSE = 0.589 and the green band of L8 producing an RMSE = 0.608, in both cases with a bias of -0.001. We found that our Bayesian inversion of PARAS resulted in LAI retrievals of slightly larger RMSE for L8. Our best result, the regularizing prior for L8, had an RMSE of 0.649, which is 0.041 units or 6.5% larger than Korhonen et al. [20] . In ad-dition, our results exhibit significantly larger biases than Korhonen et al. [20] , which originate from the fact that physically-based models attempt to generalize across sensors and viewing geometries, but are not able to account for all influences on the optical signal (e.g. atmospheric correction errors) and further need to simplify the forest radiative transfer process to facilitate computations. The comparison of our results to the empirical study of Korhonen et al. [20] may indicate the potential of reflectance model inversion for L8 images in boreal forest. On the other hand, our results for S2 were significantly less accurate than the results of Korhonen et al. [20] , both in terms of bias and RMSE. Besides the fundamental model difference between our study and Korhonen et al. [20] , who used empirical regression models to retrieve LAI for the lidar plots, they also carried out a band selection with a penalty term for multiple band predictor models. We did not aim at selecting individual bands to test their potential for LAI retrieval accuracy, as this practice seemed to contradict the physical approach, which aims at utilising the entire optical signal in an attempt to generalize the radiative transfer process in forests.
In all prior models, Landsat 8 clearly outperformed Sentinel-2 in retrieving LAI, with the bias comprising the major difference in accuracy. We attribute S2's tendency towards a higher bias in LAI retrievals to two characteristics; (i) the differences between the sensors in the visible bands, and (ii) the increased spectral resolution. The former, as discussed in detail above, is due to S2 visible reflectance being lower than L8 in our data. PARAS model has shown to have a tendency to overestimate reflectance of a given stand in boreal forests, which in turn lead to an overestimation of leaf area index [15,46] . Hence, a lower reflectance in the visible bands of S2 lead to an increased absolute bias in the LAI retrievals. In the Bayesian inversion approach, all model parameters are retrieved simultaneously. Our results indicate that the inversion approach was able to partly compensate for the lower reflectance in S2 by reducing the leaf albedo in the respective bands of conifers, which were the dominant species. However, as the spectral bands were modeled as partly correlated, the reduction in the visible bands also affected the other bands, which may have inhibited the compensation for the difference in visible bands of S2.
Modeling the spectral bands as partly correlated was necessary, as exploratory analysis of modeling the spectral bands as independent indicated an alternative solution which mistook deciduous foliage with higher albedo for understory vegetation, thus leading to substantial errors in the LAI retrievals. The second potential reason for the increased bias is the higher spectral resolution of S2. While a higher spectral resolution would be expected to yield higher accuracy of inversion results [37] , we observed the opposite to be the case. This issue may be unique to our dataset and modeling approach. The increased spectral resolution of S2 is exclusively in the red edge region. This spectral region is characterized by a steep rise in leaf albedo, from levels of about 0.1 to about 0.8 on a wavelength range of about 50 to 80 nm, constituting the steepest slope in leaf spectra in the optical region. Due to this steepness, any error, potentially originating from the input spectral data, or the sensor spectral response functions, may have had an exceedingly large influence on the results relative to other bands which are located in relatively flat spectral regions.
In addition, we transformed the spectral data to the sensorspecific bands before applying the PARAS model, which may have had a strong influence due to its nonlinearity. This a priori transformation was unavoidable for technical reasons, and was expected to not result in significant errors. However, a shift in the red edge infliction point by a few nanometers may already cause significant changes in the spectra, which were impossible to model in our approach. We carried out an exploratory analysis, where we adapted the visible bands of S2 to the levels of L8 using the linear regression equations in Fig. 1 as well as chose a subset of the bands of S2 that correspond to L8. We found that adapting the visible bands of S2 lead to a large bias reduction, and excluding the red edge bands further slightly reduced the bias, however, at best the bias of S2 got close to that of L8, but never less.
Based on our results, it appears that the saturation effect in relation to LAI is a dominant factor that is to be taken into account when inverting a forest reflectance model. Our recommendation is to regularize LAI when inverting a forest reflectance model, as our results show. The regularizing prior retrieved LAI best for L8, and only slightly worse in S2. Compared to using a realistic prior, which resulted in very similar accuracy for LAI retrieval, the regularizing prior seems to be favorable from a theoretical point of view. When using a regularizing prior, it is not necessary to obtain information on a realistic LAI distribution in the study area, but rather it is possible to define a suitable prior distribution by examining the forest radiation regime on a more general level, that is, based on radiative transfer and its characteristics in forest reflectance models. This physical foundation offers theoretically higher potential for generalization because it is only necessary to understand when the optical signal saturates, rather than having to make assumptions specific to a study site. However, previous work with a Bayesian reflectance model inversion by Varvia et al. [46] retrieved LAI using a uniform prior using EO1-Hyperion hyperspectral data. Hyperion data has been found to produce less biased estimates of leaf biochemical traits than Landsat 8 in a Bayesian inversion approach using PROSPECT model at leaf level [37] . Shiklomanov et al. [37] further reported a generally increasing bias with decreasing spectral resolution. This may indicate that when using data with low spectral resolution, the choice of priors becomes increasingly important.
The reported values for LAI retrievals follow the definition for the true LAI, i.e. the unit hemisurface area of all leaves per unit ground area. In more recent publications, the LAI has been termed plant area index, which includes the contribution of other plant elements such as woody parts. We acknowledge the importance of accounting for the contribution of non-leaf plant elements to the forest radiation regime, but used the terms LAI and plant area index interchangeably. The true LAI is the product of the effective LAI [3] and a clumping index. The effective LAI can be interpreted as the LAI of a stand with equivalent angular gap fractions but purely randomly distributed foliage, which is necessary for utilizing the modified Beer's law [35] in the forest radiation regime. The clumping index is the single factor describing the spatial aggregation, or non-randomness of the leaf distribution in a canopy. Hence, a retrieval of true LAI is necessarily a joint prediction of effective LAI and clumping index. The three prior models we used were for the effective LAI, as the clumping index was modeled separately, and the true LAI was derived from those two factors as L = L e / CI . Therefore, our results for LAI retrieval can only be interpreted in combination with the estimates for clumping index.
The posterior modes for both sensors show a similar relationship between the conifer share and clumping index of plots ( Fig. 5 ). The range of the clumping index stretches further in both directions for S2, which also exhibits a slightly wider range of conifer share values. Based on the prior formulation, the relationship between conifer share and clumping should be a linear mixture, with a clumping index of around 1 in pure deciduous stands, and 0.6 in pure conifer stands. However, our results show that the clumping index decreases stronger than the prior formulation, i.e. at a conifer share between 0.7 and 0.8, there is a drop in clumping index down to a posterior mode of slightly below 0.5, with the highest conifer shares in the posterior at around 0.83. Commonly, clumping correction is done only for conifer shoots, following Oker-Blom and Smolander [29] , while effort s to understand clumping of plant elements at higher hierarchical levels have been few [42] . While considering the overall modeling uncertainty of this study, the posterior marginal distribution of the clumping index of results data hints at increased clumping that ranges beyond shoot level, i.e. at clumping indices below 0.4, which is the lower end of observed shoot clumping [29,30] . While these findings must be interpreted with caution due to a lack of validation, the retrievals of true LAI can be indirectly used to validate the possibility of clumping beyond shoot level having a stronger influence on the forest radiation regime as described by the PARAS model than previously thought. Furthermore, these results are only reasonably obtainable by using a Bayesian inversion approach, or any other inversion approach that has inherent capabilities for uncertainty quantification.
These speculative findings, however, can only be interpreted while keeping the assumptions and potential limitations of the data, the Bayesian inversion approach, the prior model, and the PARAS reflectance model in mind. The hierarchical modeling approach, i.e. the spectral properties of leaves and understory being modeled across all stands, did not allow for variation in optical properties between stands. This modeling decision was made due to both computational and modeling reasons. Modeling optical properties at stand level would have, due to the large number of plots, required enormous amounts of computer memory that were not available. Furthermore, modeling optical properties at stand level would have potentially lead to parameter identification issues, which we encountered during an exploratory test where we modeled all parameters at stand level. The parameter identification problems were particularly prominent in the distinction between tree species, leading to confusion between coniferous and deciduous dominated stands, and subsequently to problems in clumping index and LAI retrievals. As we were particularly interested in species distinction in relation to clumping, we chose to model optical properties over the whole study area. The understory optical properties were modeled accordingly, leading to the inversion model being unable to distinguish differences in understory types. This, similarly, was necessary to avoid parameter identification problems as a consequence of the low spectral resolution of the sensors. However, this modeling choice was practical for our study, and we point out the need for further investigation on a robust hierarchical modeling framework for reflectance model inversion that utilizes synergies between individual stands and while efficiently accounting for and enabling the analysis of variability between stands.
Building on the work of Varvia et al. [45,46] , we applied the Bayesian inversion approach to multispectral data. However, adapting the approach developed for hyperspectral data required us to change their prior model. Compared to Varvia et al. [45] , we defined priors for optical properties that were less correlated, which is reasonable given the further distance in wavelength between multispectral and hyperspectral data. We further enhanced the model to take into account tree species shares at stand level, to allow a more in-depth analysis of which combinations of forest stand parameters likely constitute the remotely sensed signal. On some parameters such as the clumping index, we imposed less constraints on the prior, i.e. by truncating the distribution at values 0.1 and 1.1. Due to the lower spectral resolution, we modeled the uncertainty around the observed bottom-of-atmosphere reflectance values with an interval of 20% to account for uncertainty introduced by noise and errors from the imaging process, preprocessing of the raw sensor data, the atmospheric correction, and the PARAS reflectance model itself. We tested uncertainty parameterizations for the BoA reflectance between 5% and 25% in steps of 5%, with 20% showing the best results in terms of both numerical stability and accuracy of LAI retrievals.
A current limitation of the Bayesian approach, particularly the MCMC method we used, the no-U-turn sampler, concerns the sampling of multimodal distributions. There are MCMC methods that are able to better sample from multimodal distributions, such as inter-chain adaptation [5] or sequential Monte Carlo [7] , however, they were not implemented in the software in a way usable for our study at the time we carried out the computations. The focus of this study was, in terms of software, to build a Bayesian reflectance model inversion approach that relies on readily avail-able software to illustrate the potential for practical applications in remote sensing of vegetation. We thus refrained from implementing custom methods other than implementing the PARAS model in the Bayesian framework. The problem of sampling from multimodal distributions with MCMC was avoided in this study by using initially loose boundaries of prior distributions, and empirically analyzing the occuring multiple modes, and applying domain knowledge to identify illogical solutions to the reflectance model inversion by setting boundaries to the priors, i.e., truncating prior distributions at values that exceed an ecologically meaningful range. We expect that, with the rapid development in Bayesian computational methods and their increasing availability to scientific communities beyond mathematics and computer science, the problem of multimodality will hopefully become less prominent in the future.
Finally, there are limitations to the results presented here that originate from the PARAS model itself. The LAI retrieval bias towards overestimation as well as the mostly lower leaf and understory reflectance factors indicate that PARAS has a tendency to overestimate reflectance of a given stand. There are a number of reasons why PARAS may overestimate forest reflectance, which come from the still poor understanding of some parameters in the model, and the model formulation itself. PARAS employs a parameter for scattering uncertainty; Q ( Eq. (5) ), which is based on a semi-physical relationship between scattering in forward and backward direction of a complex medium parameterized by empirical ray tracing model in a Scots pine stand [28] . While parameterizing the canopy bidirectional scattering phase function using Q is currently the most reasonable assumption, this parameter would greatly benefit from further study. Next, it is common to correct coniferous forest LAI for clumping at shoot level, however, clumping at higher levels may play an important role, but is so far poorly studied [42] . Few studies [e.g. 24,52 ] were looking into clumping based on terrestrial laser scanning, we hope that these effort s will continue and become more extensive. Clumping, in PARAS model and, more generally, in forest radiation regime models, plays a crucial role as it governs the relationship between the true LAI and angular gap fractions, which would be the single parameter, despite the leaf angle distribution, that is able to compensate for the inconsistent estimation accuracy between LAI and canopy gaps we found in this study. Finally, though to a lesser extent, PARAS does not account for multiple interactions of photons between the canopy and understory. Multiple interactions have been studied in PARAS and found to contribute significantly to stand albedo in conditions where the understory is covered in snow, but otherwise multiple interactions may not play a strong role in forest reflectance [26] .
Besides these issues with the PARAS model, we expect, though to a lesser extent, that the assumptions around using Beer's law and the G-function for angular leaf projection area ratio may have contributed to the inconsistencies in the angular gap estimates. Beer's law, together with the G-function, which incorporates the leaf angle distributions, was used to estimate angular gap fractions based on LAI. Here again, clumping plays a significant role, as it increases or decreases the gap fraction at a given LAI and leaf angle distribution. The G-function we used for conifer trees were based on Stenberg [39] , who derived a version of the G-function for non-flat leaves. While using these shapes for conifer needles seems to be more reasonable, the application of this G-function strongly affects the common assumptions of a spherical leaf angle distribution, i.e. that G is not constant for any angle. Furthermore, the conifer G-function of Stenberg [39] has, to our knowledge, not been applied in reflectance model inversion studies. Lastly, the leaf angle distributions we used were planophile for deciduous trees, and spherical for conifers. While the planophile distribution was clearly recommended by Pisek et al. [31] for broadleaved trees in boreal forests, we did not find similarly detailed studies for coniferous trees in boreal forest.
The above limitations are mostly pointing at inconsistencies between the often necessary assumptions in forest reflectance modeling and the observations in the field. However, despite all these limitations we have shown in this study that inverting the PARAS reflectance model using multispectral satellite data is not only possible, but leads to a reasonable accuracy in the retrieval of leaf area index while simultaneously quantifying the uncertainty related to the retrieval. The latter quality, the uncertainty quantification, is in our opinion crucial, particularly in the presence of model bias and the optical saturation effect, as inversion methods that do not provide information on uncertainty may lead to false intuitions about the underlying retrieval accuracy, which further leads to false expectations of potential users of inversion methods, and subsequently may lead to a biased interpretation of the utility of forest reflectance model inversion in forest biophysical parameter retrieval.

Conclusion
In this study, we have developed a Bayesian reflectance model inversion approach for multispectral satellite images and applied it to two medium resolution satellite images from Sentinel-2 and Landsat 8. We found that, in multispectral images, the inversion of a reflectance model greatly benefits from regularization in the prior model, and that the inversion approach was able to produce reasonable estimates of leaf area index without relying on any further information from the study site other than the satellite images. We further discussed the limitations of reflectance modeling with the example of PARAS model, the Bayesian inversion approach and their respective limitations, and pointed out existing gaps in our understanding of the forest radiation regime as expressed through common modeling assumptions.