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

ContentslistsavailableatScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer

journalhomepage:www.elsevier.com/locate/jqsrt

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

Daniel Schraik

^a,∗

, Petri Varvia

^b

, Lauri Korhonen

^c

, Miina Rautiainen

^a,d

aDepartment of Built Environment, School of Engineering, Aalto University, Finland

bFaculty of Information Technology and Communication Sciences, Tampere University, Finland

cSchool of Forest Sciences, University of Eastern Finland, Finland

dDepartment of Electronics and Nanoengineering, School of Electrical Engineering, Aalto University, Finland

a rt i c l e i n f o

Article history:

Received 11 March 2019 Revised 13 May 2019 Accepted 13 May 2019 Available online 14 May 2019 Keywords:

Forest reflectance Bayesian inversion PARAS

Sentinel-2 Landsat 8 Leaf area index Clumping Spectral invariants Recollision probability

a b s t r a c t

Theinversionofreflectancemodelsisageneralizabletooltoobtainestimatesonforestbiophysicalpa- rameters,suchasleafareaindex,with theoreticallylittleinformation needfromastudyarea,instead relyingontheknowledgeaboutphysicalprocessesintheforestradiationregime.Theuseofpriorinfor- mationcangreatlyimprovethereflectancemodelinversion,however,theliteraturedoesnotyetprovide muchinformationontheselectionofpriorsandtheirinfluenceontheinversionresults.Inthisstudy,we usedaBayesianapproachtoinvertthePARASforestreflectancemodelandretrieveleafareaindexfrom Sentinel-2MSIandLandsat8OLImultispectralsatelliteimages.ThePARASmodelisbasedonthethe- oryofspectralinvariants,whichdescribestheinfluenceofwavelength-independentparametersonforest radiativetransfer.TheBayesianinversionapproachishighlyflexible,providesuncertaintyquantification, andenablestheexplicitincorporationofpriorknowledgeintotheinversionprocess.Wefoundthatthe choiceofpriorinformationiscrucialininvertingaforestreflectance modeltopredictleafareaindex.

Regularizingandinformativepriorsforleafareaindexstronglyimprovedthepredictions,relativetoan uninformativeprior,inthattheycounteractedthesaturationeffectoftheopticalsignaloccuringathigh valuesforleafareaindex.ThepredictionsofleafareaindexweremoreaccurateforLandsat8thanfor Sentinel-2,duetopotentialinconsistenciesinthevisiblebandsofSentinel-2inourdata,andthehigher spectralresolution.

© 2019TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Forest reflectance models aim atproviding 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 havebeensuggestedinthe literature[e.g.19],andmostmethods producepoint estimatesasoutputoftheforest reflectancemodel inversion. These point estimates are obtained by applying one of several radiative transfer inversion methods, that are, look-up tables,optimizationmethods,andneural networks,amongothers.

Theseinversion methodshavebeenproven toprovide reasonable retrieval accuracy of forest biophysical variables, given that the

∗ Corresponding author.

E-mail address: daniel.schraik@aalto.fi(D. Schraik).

constructionofthelook-uptable,orthetrainingdataforaneural network, cover the ecologically reasonable parameter range of radiative transfer model inputs [19].However, the accuracy with whichmanyinversionmethodspredictforestbiophysicalparame- tersstronglydependsonthewaytheyareparameterized,orwhat initial guess one chooses. Optimization algorithms in particular maybe challenging dueto their need for an initial value guess, whichcanleadtoasignificanterror[45].Recentadvancesincom- putationsandalgorithms relatedtoBayesianinferencehavemade 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, thoseunknownsare theforeststructuralandspectralparameters, conditional on observed bottom-of-atmosphere (BoA)reflectance.

The Bayesian framework features, among others, two properties thatmakeitparticularlysuitableforreflectance modelinversions, (1)theinversionresultbeingaprobabilitydistributionratherthan https://doi.org/10.1016/j.jqsrt.2019.05.013

0022-4073/© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )

a point estimate, which provides inherent uncertainty quantifi- cation, and (2) the possibility to explicitly incorporate existing knowledgeonparametersintothemodelitself.Afewstudieshave appliedaBayesianinversionapproachtoforestreflectancemodels [45,46,51]ortheleafopticalmodelPROSPECT[37]withpromising results,particularly interms of uncertainty quantification.Varvia etal.[45] haslaid out the groundwork fora Bayesian inversion approachusingthePARASforestreflectancemodel[34].

The PARAS model is based on the theory of spectral invari- ants,whichdescribestheinfluenceofcanopystructureonthefor- estradiationregime.Themodelparameterizestheforest radiation regime using only few spectral and structuralparameters, which makes it particularly suitable for a Bayesian inversion. Further- more,PARAShasalreadyshownpromisingresultsintheinversion ofhyperspectral datatoretrieveleafareaindex, andleaf andun- derstoryopticalproperties[45,46].

Physically-based reflectance model inversion usually sets out to retrieveforest biophysical parameters, such asleaf area index or other canopy structural or biochemical parameters, from ob- served airborneor spaceborne imagery. Thisbiophysical retrieval approachis only in part competingwith 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 empiri- cal approach ismore suitable forregional applications withhigh predictionaccuracyrequirements[38].Akeyproperty ofphysical modelsistheirattempttobesensorindependent,whichisbecom- ing more important with every new satellite mission to be able tomaintainaconsistenttheoreticalfoundationforcombiningdata fromdifferentsensors.

In optical remote sensing of vegetation, two recent satellite missions have attracted a lot of attention in science and indus- trycommunities;the highresolution optical satellite sensorson- board Sentinel-2a and -2b (Multi Spectral Instrument, MSI), and onboardLandsat 8 (Operational LandImager, OLI). Hereafter, we usethe namesofthesatellite platformsandtheir respectivesen- sors synonymously, abbreviated as S2 and L8. These two multi- spectralsensors share bandsofsimilar spectral sensitivity,which appeartogenerallyagreewell [e.g.[4,22,25]], despitesome stud- iesreportingdifferencesparticularlyintheredband,whereS2ap- pearstoexhibitlowerreflectanceinsomecases[20,44,47].S2fur- therfeatures three bandsin the red edge spectral region, which hasbeen shownto slightly improveretrievals of forest biophysi- calparameters[20].Besidesthenoveltyoftheadditionalrededge bandsofS2,bothsensorshaveanunprecedentedspatialresolution forfreely available satellite imagery; 10/20mresolution forvege- tationbands of S2, and 30m resolution for L8. From a Bayesian reflectancemodelinversionperspective,thesetwosensorsare in- terestingbecauseprevious studiesusedPARASmodelonlyonhy- perspectraldata[46].The effectof spectral resolutionwasprevi- ouslystudiedby Shiklomanovetal.[37]atleaflevel,wherethey pointedout that inversion uncertainty increases withdecreasing spectralresolution. Theyalsonotedthat L8in particularhas,de- spiteonlysubtledifferencesinspectralresponse,performedbetter inthemodelinversionthanitspredecessorsLandsat5andLandsat 7,potentially hintingthatthemultispectral resolutionofL8tobe inaparticularlyinterestingrange.

Despitemanyeffortsforsensor-intercomparisoninforestappli- cationsingeneral[e.g. 1,12],studies thatanalyze therelative po- tentialofS2 andL8for forest biophysicalparameter retrieval re- mainfew,andarelimitedto empiricalmodeling approaches[e.g.

20].Toourknowledge,thisisthefirststudytoapplyaphysically- basedforest biophysicalparameterretrievalusingS2andL8data.

BesidestheattempttocompareS2andL8intheirpotentialforfor- estreflectancemodelinversion,wefurthersetout(1)todetermine whetherthisinversionapproachisfeasibleatthislevelofspectral

resolution,(2)toestimate theinversionqualitythroughtheaccu- racyofleafarea indexretrievals,and(3)toanalyzetheinfluence ofthechoiceofpriordistributionontheinversionaccuracy.

2. Materialandmethods 2.1. Studyarea,fieldandlidarplots

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 sylvestrisL.),andSilverandDownybirch(BetulapendulaRothand Betula pubescens Ehrh.). The area is characterized by diverse for- est structure and density, with average stand size of about 2ha.

Ourreferencedataconsistedof746lidar plotsforwhich theleaf area index(L)waspredictedusingairbornelidar data.TheLiDAR datawerecollectedon4September,2014usingaLeicaALS70-HA laserscanningsystem.Ithasawavelengthof1064nmandabeam divergence of 0.22 mrad. The flight altitude was 2000m above groundlevel,andthescananglewasconstrainedto ± 15degrees.

The nominal pulsedensity ofthe lidar data was0.75 pulses per squaremeter,withafootprintsizeofabout44cm.Therewere19 plots withfield-measuredLwithin the scanned area which were usedtotrainregressionmodelsforpredictingLforthelidarplots.

L was retrieved from the fraction of all echos over 1.3m above ground, andcanopycoverfromthe fractionofallfirst andsingle echosover1.3maboveground,basedonKorhonenetal.[21].The lidarplotshadasquareshapewithedgelengthof60m,andwere placedovera setofnine(S2)orfour(L8)pixels.Dueto geoloca- tiondifferencesoftheimages,theplotcenterpointswerelocated around 16m away fromeach other.More details are available in Korhonenetal.[20],wherethedatasetwasfirstdescribedandis availablefromtheirsupplementalmaterials.

2.2. Satelliteimages

We used multispectral satellite images from Sentinel-2a MSI (S2)andLandsat8OLI(L8)asprovidedinthedatasetofKorhonen etal.[20].TheS2andL8imageswereacquiredon17August2015 09:50amUTC(absoluteorbitnumber000791,tilenumber35VNK, solar zenith angle 49.3^◦), and on 22 August 2015 09:22am UTC (path187,row16,solarzenithangle51.7^◦),respectively,temporally closetothefieldcampaign.

The S2 image was obtained from Sentinel Scientific Data Hub [10] as Level-1C product (geometrically corrected top-of- atmospherereflectance).Anatmosphericcorrection wasdoneus- ingtheSen2cormoduleversion2.2.1withintheSentinel2Toolbox (S2TBX),SentinelApplicationPlatform(SNAP)version4.0.2.Inthis study, spectral bands with 10m and 20m resolution were used.

S2’s band 8 was omitted from the analysis due to potential in- consistencieswithL8,asrecommendedby Lietal.[22].Geomet- rically corrected Landsat 8 Surface (bottom-of-atmosphere, BoA) ReflectanceClimateData Records[27]were downloadedfromthe UnitedStatesGeologicalSurvey[43].Thesquare-shapedlidarplots containedtheaverageofnine(S2)andfour(L8)pixels.

2.3. PARASmodelandparameterization

TheforestreflectancemodelPARAS[34]isbasedontheconcept ofspectralinvariants,specificallythephotonrecollisionprobability, andrelies onlyon a handfulofforest structural andspectral pa- rameters.ThemodelsimplicitymakesPARASsuitableforinversion withMarkovChainMonteCarlo,astheill-posednessofaradiative transfermodelisexpectedtoincreasewiththenumberofparam- eters,asparameter identificationbecomesmoreproblematic[e.g.

32].PARASmodelsforestreflectanceas

r

( θ

i,

θ

v,

λ )

^=rG

( θ

i,

θ

v,

λ )

^+rC

( θ

i,

θ

v,

λ )

^{, (1)}

orinother words,theforest standbidirectionalreflectancefactor (r) consistsof a reflectance contribution fromthe groundor un- derstory (r_G)andthetreecanopy (r_C). Bothcontributionsdepend ontheilluminationandviewingzenithangles(

θ

iand

θ

v,respec- tively),andthewavelength(

λ

^).

The understorycontributiondescribesthe proportionofsunlit groundareathatisvisibletothesensor:

r_G

( θ

i,

θ

v,

λ )

=T

( θ

i

)

^T

( θ

v

) ρ

G

( λ )

, (2) with T, canopy transmission (also called canopy gap fraction), forming the bidirectional gap probability T(

θ

i)T(

θ

^{v) in illumina-}

tion and viewing direction, and the understory hemispherical- directionalreflectancefactor(

ρ

G).

Thetreecanopycontributiontothereflectedsignalismodeled as the fraction of intercepted sunlight that is scattered into the viewingdirection:

rC

( θ

i,

θ

v,

λ )

=i0

( θ

i

) ω

C

( λ )

^Q

( λ )

, (3) with the canopy interceptance i₀(

θ

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- metricseries

ω

C=

ω

L(¹−p)+

ω

²Lp(¹−p)+...+

ω

ⁿL⁺¹pⁿ(¹−p)=

_∞

k=1

ω

^kL⁺¹p^k(¹−p),wherekdenotesthescatteringorder(i.e.k= 1isfirst-orderscattering,orphotonsthat escapethecanopyafter oneinteraction).Thisgeometricseriesconvergesto

ω

^C

( λ )

=

ω

^L

( λ )(

¹−p

)

_1−p¹

ω

^L

( λ )

⁽⁴⁾

withthecommonratio

ω

Lp,asthephotonhasadditionalinterac- tionswithcanopyelementsandisscattered,andthefirsttermof theseriesis

ω

L(¹−p),which,fromacausalperspective,areactu- allythefirstandlast partsofthesequence.Thatis,thephotonis scatteredfromitsinitialinteractionwithinthecanopy(

ω

L)andis ultimatelyescapingthecanopy(1−p).

Qisthedirectionalparameterdescribingthefractionofcanopy- scattered radiation that escapes in an upward direction, i.e. to- wards the sensor. This backscatterfraction wasparameterised as proposedbyMõttusandStenberg[28]:

Q=1 2

1+q1−p

ω

^L

( λ )

1−pq

ω

^L

( λ )

. (5)

Thisparameterisation ofQ incorporatesan additional spectralin- variantparameter,q,whichwasempiricallydeterminedbyMõttus and Stenberg [28] asq=exp(−0.1684Le)^{. The effective leaf area} indexLe inthiscontext isrelatedtothe (true)leafarea indexas Le=

β

^{Lvia theclumpingindex}

β

^{.We usethetermeffectiveleaf}

areaindexasdefinedinBlacketal.[3].

Thephotonrecollisionprobabilitypcontainsthemostinforma- tion onforest structure inPARAS model,anditdependson both clumpingandtheangulardependencyofcanopygapsovertheen- tirehemisphere.WeparameterizedpfollowingStenberg[40]: p=1−iD

β

Le , (6)

with the diffuseinterceptance i_D,the clumping index

β

^andthe

effectiveleaf area indexL_e.The diffuseinterceptanceinPARAS is describedas

i_D=1−2 _π/2

0

T

( θ )

^sin

θ

^cos

θ

^d

θ

. (7)

Theangularrelationship betweenLe andTwasmodeledusing themodifiedBeer’slawwithaconstantclumpingindex:

T

( θ )

=exp

−G

( θ )

^Le

cos

θ

, (8)

withtheprojectionofunit leafarea intoa zenithangleG(

θ

^).Pa-

rameterizing G requires estimates of the leaf angle distribution, whichweincorporatedinourinversionmodel.Gdependsonleaf angle(

θ

L)as[49]

G

( θ )

= _π/2

0

A

( θ

^,

θ

^L

)

^f

( θ

^L

)

^d

θ

^{L. (9)}

Ainthisequationdescribestheprojectionarearatioofleaveswith zenith angle

θ

L into the direction

θ

^{, and} incorporates assump- tions on the shape of leaves,and f is the frequency distribution ofleafangles.Wemodeled Aseparatelyfordeciduousandconifer trees.Fordeciduousleaves,weusedtheprojectionfunctionforflat leaves[48]

A

( θ

,

θ

L

)

=

cos

θ

^cos

θ

L,

|

^cot

θ

^cot

θ

L

|

>1 cos

θ

^cos

θ

L 1+_π²

(

^tan

−

)

,otherwise, (10) with=cos⁻¹(^cot

θ

^cot

θ

L)^{.Forconifer leaves,weusedthefunc-} tionfornon-flatleavesproposedbyStenberg[39]:

A

( θ

^,

θ

^L

)

⁼

π

²²

_π

0

1−

(

^sin

θ

^sin

θ

^L+cos

θ

^cos

θ

^Lcos

τ )

^2d

τ

^{. (11)}

Wemodeledfalsoseparatelyfordeciduousandconiferleaves.De- ciduoustreesinFinlandarepredominantlybirchspecies,forwhich datawasavailablefromPiseketal.[31],whoalsoclearlyrecom- mended the use of planophile leaf angle distributions for boreal broadleavedtreespecies.Therefore,weusedatwo-parameterBeta distribution[13]tomodeltheleafangledistributionofdeciduous leavesfromtheexistingdataas:

f

( θ

^L

)

= 1 B

( μ

^,

ν )

1−

θ

^L

90

μ−1

θ

^L

90

ν−1

, (12)

withtheBetafunctionB,andtwoparameters

μ

^and

ν

^,whichde-

pendonthemeanandstandarddeviationof

θ

L.Thisdistribution wasrecommended by Wang etal. [48] for cases where data on meanandstandarddeviationisavailable.Therewasalackofdata forboreal conifer tree species leaf angle distributions, hence we assumedasphericalleafangledistribution,followingDeWit[6]

f

( θ

^L

)

=sin

( θ

^L

)

. (13)

WeapproximatedtheintegralsinEqs.(9)–(11)byquadratureusing 10 evenly spaced locations, and Eqs. (7) and (8) at five discrete angleswithassociatedweights,correspondingtothemethodused fortheLAI-2000instrument[23].

In a nutshell, the PARAS model has the unknown parameters

ρ

G,

ω

L, T(

θ

i), T(

θ

^{v), q, p.The latter two} parameters, q andp,re- quirefurtherparameterization,andinparticularprequiresacom- pleteparameterizationofforest structure,namelythroughtheef- fectiveleafareaindex,theclumpingindex,andtheleafangledis- tribution.Weassumedaplanophileleafangledistributionforde- ciduous trees [31], and a spherical distribution for conifer trees.

ThePARASmodeldoesnotinherentlymodeltreespecies,however, itis possible tomodel treespeciesthrough a mixture parameter forleaf optical properties andclumping index. We distinguished onlydeciduousandconifertreespeciestoavoidparameteridenti- ficationproblemsbetweenthetwoconiferspecies,Norwayspruce andScotspine,andrefertothespeciesparameterasconifershare, which is definedas the share ofconifer speciesin the leaf area index. Parameters that are modeled by speciesare denoted by a subscript“c” forconifer,and“d” fordeciduous.

2.4.Bayesianreflectancemodelinversion

TheforwardproblemofPARAS,whichsimulatesforeststandre- flectancefromgivenforest structuralandspectral parameters,can

bestatedasr=h(^x)+

^{.Inthiscontext, r}∈R^m×n istheobserva- tion matrix, it contains the satellite pixel measurements for the numberofplots m=746 andthe numberof sensorbands n,i.e.

n=6forL8,andn=9forS2.

Thevectorofunknownparametersx=[

ρ

G^T,λ,

ω

^TL,c,λ,

ω

^TL,d,λ,L_e,i,

β

c,i,

β

d,i,p_c,i]^T∈R^4m+3ncontainsallinputparametersofthePARAS modelina hierachical modeling approach.The optical properties

ρ

G,

ω

L,c, and

ω

L,d were modeled foreach spectral band

λ

^{of the}

respectivesensor,atthelevel oftheentirestudyarea. Theforest structuralpropertiesLe,

β

^c,

β

d,pcweremodeledforeachplotiin thestudy area. We chose this hierarchicalmodeling approach to leverageassumedsimilaritiesbetweenplots.

The following elaboration on theBayesian inversion of PARAS closelyfollows Varvia et al.[45]. Inthe Bayesian approach,both observedandunknown parameters are modeled asrandom vari- ables. We therefore denote the prior probability density by P(x), which contains available information on x before any measure- mentsweredone.Wereliedonreadilyavailableinformationfrom peer-reviewedliteraturein combinationwithsome practical con- straints to select suitable prior distributions for each parameter inx. More details on prior selection are given in Section 2.5. In Bayesianinference,thepriordensityisupdatedbytheinformation gainedfromsatellitemeasurementsusingBayestheorem

P

(

^x

|

^r

)

=P

(

^r

|

^x

)

^P

(

^x

)

P

(

^r

)

^∝P

(

^r

|

^x

)

^P

(

^x

)

, (14) whereP(x|r)istheposteriorprobabilitydensity,P(r|x)isthelikeli- hoodfunction,andP(r) istheposteriormarginaldensity.Thelat- ter termserves to normalizethe posterior, andcan be neglected inthefurtheranalysis.Theposteriordensityisthefullsolutionto theBayesianinverseproblems,i.e.itisa jointprobability density ofallforeststructuralandspectralparametersinPARAS.Thelike- lihood function containsthe information fromthe measurement, i.e.the satellite bands’observedBoA reflectance factorfortheli- darplots. The likelihoodis coincidentwiththe forwardproblem, i.e.r=h(^x)+

,where

∈Rⁿisanadditiveerrorterm,containing bothmodelerrorandmeasurementnoise. Inthecaseofadditive error,thelikelihoodfunctionbecomes

P

(

^r

|

^x

)

=P

(

^r−h

(

^x

))

, (15) withtheerrordensityP modeledasamultivariatenormaldistri- butionwithzeromean,andanuncorrelatedcovariancewithstan- darddeviationof20%oftheobservationsineachband.

2.5.Priorchoices

Tostudytheeffectofpriorchoicesontheinversionresult,we chosethreedifferentprior distributions forthe effectiveleafarea index:

P

(

^Le

)

=

1/10, 0≤Le≤10

0, otherwise , (16)

P

(

^Le

)

=

N

(

²,1

)

, 0≤Le≤10

0, otherwise ,and (17)

P

(

^Le

)

=

N

(

⁰,2

)

, 0≤Le≤10

0, otherwise , (18)

where N denotes a normal distribution, with their respective meansandstandarddeviations.

Thesethreepriorsresembledifferentmodelingapproaches.The firstapproachisuninformative(Eq.(16)),itmakesnoassumption aboutLe other thanimposingarangeofpossiblevalues.Thesec- ondpriorisinformative(Eq.(17)),itassumesthattheprobability distributionforLe inanyplotfollowsthedistributionofLe inthe

geographical region. We approximatedthis assumption by a nor- maldistribution,basedtheresultsofHärkönenetal.[14],whoes- timatedLe fromtheFinnishNationalForestInventory data2004–

2008. The third prior is regularizing (Eq.(17)), it formulates the probabilitydistributionofLeparticularlytoregularizetheinversion inawaythatcounteractsasaturationeffectoftheopticalsignalat highvaluesofL_e.Thissaturationeffectcausestherelationshipbe- tweenLeandforeststandreflectancetoflattenoutindenseforest stands,leadingtohighuncertaintyandpotentialoverestimationof L_e in the inversion(Fig.3). Hence, the third prior waschosen in anticipationandtocounteractsuchasaturationeffect.

The prior for conifer clumping index (

β

^{) is based on the as-}

sumptionthat shoot clumpingisthe main causefornon-random distribution ofcanopy elements. FollowingStenberg etal.[41], a shootclumpingindexcanbeexpressedasfourtimesthesilhouette tototalarearatio(STAR)ofashoot.Theobservedrangeforshoot STARranges fromabout 0.1to 0.22[41],resulting ina clumping indexof 0.4to 0.88.For theprior of

β

^{,we chose a meanof 0.6}

anda standard deviationof 0.2to assign the highest probability totheobservedrangeofshootSTAR,butsimultaneouslyallowthe clumpingindextotakevaluesbeyondthisrange,whichreflectson thelittle knowninfluence ofhigherlevel clumping.In deciduous trees, shoot clumping is non-existent, and the common assump- tionisthattherandomleafdistributionholds.Tobeabletoquan- tifyapotentialinfluenceonhigherlevelclumping,e.g.atcrownor standlevel,we modeledaseparate clumpingindexfordeciduous trees centeredat1, withastandard deviationof 0.2. Theclump- ingindexforbothspecieswereconstrainedwithintherange0.05 and1.1.Thedominanttreespecieswasassumedtobe coniferous, expressedthroughthepriorfortheconifersharepcwithmeanof 0.8andstandarddeviationof0.5.

Priorsforleafandunderstoryopticalpropertieswereobtained from published measurement data. We used the spectrometer measurementsfromHovietal.[17]forpine,spruce,andbirchleaf optical properties.Thisdatabase containsseasonalmeasurements of25borealtreespecies,withreflectanceandtransmittancemea- surementsreportedforseveralsamplesandstratifiedbyshootage, leaf side,and crownposition.We selected themeasurements for birchon 25August2016,andforpine andspruceon12 Septem- ber 2016,which were theclosest dates to thedayof yearofthe satellite image acquisitions. We added reflectance and transmit- tance measurements to get leaf albedo, andaveraged over shoot age,leafside,andcrownpositionaswell asoverindividualsam- ples.Finally,thespectraforpineandsprucewereaveragedtoform a single conifer leaf albedo. Understoryhemispherical-directional reflectance factor (HDRF)data came fromHoviet al.[16],where weusedthemeasurementsfromJune2010,andaveragedoverthe sub-xeric, mesic, andherb rich sitetypes. The spectral datawas availableinthewavelengthregionbetween400and2400nm.We convertedthespectratotherespectivesensorbandsofS2andL8 using the sensor-specific spectral response functions. Thesewere obtainedfromtheSentinel2DocumentLibrary[9]andBarsietal.

[2].Theresultingmeanspectraforleafandunderstorywereused asthepriormeansforthedeciduous,conifer,andunderstoryopti- calproperties.

Spectral properties were modeled as multivariate normal dis- tributions, centered attheir measured values.The covariance be- tweentheindividualbandswascrucialforforestreflectancemodel inversiontoavoidparameteridentificationproblems.Increatinga covariancematrix,wefollowedthecorrelationconceptlaidoutby Varvia etal.[45].The correlationmatrixR betweenbandscomes fromthreeweightedcomponents:

R=1

κ

all+Rpart

κ

part+I

κ

ind, (19) wheretheunitarymatrixmeansthecorrelationbetweenallbands, Rpartmeansthecorrelationbetweenselectedgroupsofbandsthat

areassumedsimilar,andtheidentitymatrixmeansarandomcom- ponentforeachindividual band.Theabove matriceshavedimen- sions n_bands×n_bands,i.e. 6×6for L8and 9×9for S2.The factors

κ

^{are the weighting} coefficients, and we used

κ

all=0.1,

κ

part= 0.2,

κ

ind.=0.7. To obtain Rpart, the sensor bands were grouped, withthevisiblebands,andtherededge1bandofS2formingthe first group,thenearinfrared(NIR)band beingits owngroup,to- gether withred edge2and3bandsofS2, andtheshortwave in- frared (SWIR) bands1 and2 forming the third group.Therefore,

Rpart=

⎡

⎢ ⎣

41×4 0 0

0 1

3×3 0

0 0 1

2×2

⎤

⎥ ⎦

⁽²⁰⁾

forS2,and,forL8,analogouslythethreepartialmatriceshavedi- mension3×3, 1,and2×2.Thecovariancematrix isobtainedby =SRS^T,whereRisthecorrelationmatrix(Eq.(19)),andSisthe diagonalmatrixofstandard deviationsofindividualbands,which weparameterizedwith10%ofthemeasuredspectral dataofeach band.

2.6. Computationsandaccuracyassessment

Foreachofthethreemodels,i.e.thedifferentpriorsforLe,and persensor,wesampled5000samplesineachof16parallelchains, withaburn-inphaseof5000samples.Thecomputationswerecar- riedoutinPythonusingthePyMC3package[36].

Weusedasetofdifferentaccuracymetricstoassesstheaccu- racy of parameter retrievals fromthe inversion. The full solution totheinverseproblemistheposteriordistribution,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 posteriormarginal distribu- tionthatcontains95%oftheprobabilitymassforagivenparame- ter.Theposteriormarginaldistributionisapproximatedbytheset ofsamplesonlyforasingleparameter,ignoringthejointprobabil- itydistribution withother parameters. Thereis apotentially infi- nite numberofintervalsthatcontain 95%oftheprobabilitymass, hence we use the commonapproach to select the narrowestin- terval,calledthehighestposteriordensity(HPD)interval.TheHPD hasastraight-forwardinterpretation,itdescribestheintervalthat contains the true value with a probability of 95%. The accuracy metric associatedwiththe HPDis thepercentageof plotswhose observedvalueslaywithintheHPDinterval.

Toallowcomparisonwithconventional,non-Bayesianaccuracy metrics based on point estimates, we used the posterior mode, i.e.the mostlikelyvalue ofthe posterior.The conventionalaccu- racymetricswereportaretheroot-meansquareerror(RMSE),the meanbias,andthebiascorrectedRMSE(cRMSE),calculatedasfol- lows:

RMSE=

n i=0

(

^yi−yˆ_i

)

²

n (21)

Bias= n

i=0

yi−yˆi

n (22)

cRMSE=

(

^yi−

(

^yˆi+bias

))

²

n , (23)

whereyandyˆaretheobservedandpredictedvariable,andnisthe numberofplots.The Pythoncodeusedforthisstudyisavailable fromthesupplementaryfile(S1).

3. Resultsanddiscussion

Wepresentourfindings intwoparts.Thefirstisadescriptive analysisoftheavailabledata,whichisessentialintheinterpreta- tionofthelatterpart,whichpresentstheresultsofourinversion andadiscussionontheinfluenceofpriorinformationontheinver- sionresults.Followingthepresentationanddiscussionofconcrete resultsisa generaldiscussion onthe potential andlimitationsof ourinversionapproach.

3.1. Descriptivesatelliteandplotdataanalysis

The satellite observations for the respective lidar plots of S2 andL8 were mostlyin agreement (Fig. 1). However, the BoA re- flectance values in the visible bands differ strongly between the sensors.TheblueandgreenbandsofS2showtheirstrongestdif- ferencesin the darker plots, withthe difference decreasing with increasing reflectance. In the red band, the reflectance values of S2areoverallshiftedbyabout0.01units,theregressionlineindi- catescloseragreementatareflectanceofabout0.1, whichiswell beyondthehighestreflectance observationinthestudyarea. Ko- rhonenetal.[20]alreadypointedoutthesereflectancedifferences betweensensors,andnotedthatthedifferenceslikelydonotstem from the atmospheric correction procedure. Furthermore, Korho- nen etal. [20] pointed out that two other studies [44,47] found similarreflectancedifferencesandspeculatingatsensordifferences asthemainreason.Laterstudies[e.g.22,50]foundthedifferences betweenthesensorstobemuchsmaller.TheESASentinel-2Team hasidentifiedafew productanomaliesduringthecommissioning periodfromJune to December 2015,andthe dataproducts were reprocessedtwice.Thefirstreprocessingwasdoneusingprocess- ingbaseline02.01,and,aftersomeproductanomalieswerefound, a second reprocessing was done using baseline 02.04, as docu- mentedintheSentinel-2MissionStatusReport44,July2016,avail- ablefromEuropean SpaceAgency [11].We comparedthe S2 L1C productsproducedbybothprocessingbaselines,andfoundnodif- ferenceintheimages.Sincethesensordifferencesinvisiblebands were only found in studies using early Sentinel-2a images (until May2016),butnotinlaterstudies,wheretheagreementbetween S2andL8seemsmuchbetter,wesuspectthatpre-processingofS2 imagesmay havebeendifferent duringits commissioning phase, orother technicalreasons mayhavebeen involved.TheBayesian reflectancemodelinversion approachwasabletocompensate for the sensor differences in the visible bands to some extent. In theposterior distribution, theleaf albedo wassignificantly lower for coniferous trees in the visible bands of S2 (Fig. 2). Conifer- ousstandstypically constitutethedarkestpixelsduetoincreased clumping[34],lowerreflectanceobservationsmayleadtoanover- estimatedconifershareintheinversion,andtoalowerleafalbedo ofconifersasthedominantspecies.

The relationship between the leaf area index, andreflectance showsa typicaldecreasingreflectancewithincreasingLAI, witha saturationeffectatLAIgreater thanabout3 to4(Fig.3). Similar totheabovediscussion,thereflectancediffersbetweenthesensors inthevisiblebands.TheforestreflectanceofstandswithsmallLAI consistsmostlyof thereflectedsignal fromthebrightunderstory [8,33].Hence, inaccordancewiththe above discussionon sensor differences,thereislittledifferenceintheblueandgreenbandfor sparsestands,i.e.plotswithLAIclosetozero.AstheLAIincreases, the reflectance of the visible bands decreases and levels out at higherlevels forL8than S2,a differencewhich isnot presentin NIRandSWIRbands.

3.2.Reflectancemodelinversionresults

Theuseof priorinformationin reflectancemodelinversion in general has been widely accepted [e.g. 19]. However, the use of

Fig. 1. Bottom-of-atmosphere reflectance of the common bands of Sentinel-2a MSI and Landsat 8 OLI. The solid line has unit inclination, and the dashed line shows the fitted linear regression line, whose formula and coefficient of determination are annotated in the respective subfigures.

Fig. 2. Leaf and understory optical properties resulting from the reflectance model inversion using the regularizing prior. S2 is in blue, L8 in red. Pale color shows the prior for optical properties, and the darker colors show the posterior marginal distributions.

prior information is usually limited to defining the range of re- flectancemodelparametersthatisecologicallymeaningful,orus- ingrandomdrawsfromafielddataset.IntermsofBayesianinfer- ence,settingarangeisequivalenttoassigningauniformdistribu- tiontoreflectancemodelparameters.Wesuspectthatthismaynot be the best choice, particularly for estimating parameters which exhibit a saturation effect with the remotely sensed signal, such asLAI,becausesuch nonlinearitiesinthemodelaffecttheproper choice of an objective, uninformativeprior [e.g.18]. We hypoth-

esizedthat usingnon-uniform priorsimprovespredictions ofop- ticallysaturatingparameters, that are,parameters whichatsome rangecannotbereliablypredictedfromtheopticalsignal,orwhich can onlybe predictedwithimmense uncertainty.Totest thishy- pothesis,weemployedthreereflectancemodelinversionstrategies that rely on different prior parameterizations for LAI. The three strategies,orpriormodels,were(i)auniformprioraspercommon assumption,(ii)arealistic(informative)priorwhichapproximated theLAIofFinnishforests[basedonNationalForestInventoryData,

Fig. 3. The relationship between bottom-of-atmosphere reflectance of S2 (blue circles) and L8 (red diagonal crosses) with leaf area index in different spectral bands of the satellite sensors.

14],and(iii)aregularizingpriorwhichwasaimedatcounteracting the saturation effectof theoptical signal at highLAIby favoring smallvaluesforLAI.

Overall,boththeinformativeandregularizingpriormodelsper- formed satisfactory, while the uniform prior was heavily biased (Fig.4).ForL8, theregularizingprior resultedin thebestLAIre- trievals (bias = 0.03, RMSE=0.65), comparedto theinformative prior (bias=−0.19,RMSE=0.67).Whilebothpriorsappearvery similarintheir predictionaccuracyofstandswithLAIlargerthan about2,thereis anotable differenceinsparse standswithsmall LAIbetween0and2.Inthesesparsestands,theinformativeprior introduced avisiblebias,causingoverestimatedLAIsespeciallyin therangeofLAIvaluesbetween0and2.Theregularizingprioras-

signsthehighestprobability tolow LAIvalues,hencethepredic- tionaccuracyishigherforlowLAIvalues,whilenotcompromising predictionsofstandswithlargeLAI.

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 isno visibledifference forthetwopriorsforLAIvaluesbetween0and2,aswasobserved inL8retrievals.Thismaybeduetothelowerreflectancevaluesof S2inthevisiblebands,whichinturncauseoverestimationofLAI intheinversion andhavemoreinfluenceon theresultsthan the differenceinLAIpriors.Theconfidenceintervalmetric,theratioof plotswithobserved LAIwithin the HPD,exhibits thesame rank- ingofpriorsasthepointestimate-based metricsinboth sensors.

Fig. 4. Observed and retrieved leaf area index. The point estimates of LAI retrievals are shown together with their 95% credible interval. Black lines (intervals) and diagonal crosses (point estimates) mean the credible interval overlaps with the observed value, red lines (intervals) and circles (point estimates) indicate no overlap.

InS2,theinformativeandregularizingpriorsrespectivelycaptured 66%and63% ofthe plotscorrectly. 87%and85% were thecorre- spondingfiguresfortheregularizingandinformativepriorsforL8, respectively. When considering uncertainty, as the share ofplots withintheHPDsshowed,thedifference betweentheregularizing andinformative priorsappears less pronounced than with point estimatemetrics.

For both sensors, the uninformative priors resulted in signif- icantly worse LAI retrievals (L8 bias = −2.67, RMSE = 2.98; S2 bias = −3.89, RMSE = 4.39). The uninformative case showed a stronginfluenceoftheopticalsignal saturationeffect,asthebias increaseslinearlywithincreasingLAI(seeFig.4).Boththeregular- izingandinformativepriorsoutperformedtheuninformativeprior byone ortwoorders ofmagnitudeinbias.The confidenceinter- valsforLAIin theuniformprior modelare relativelywide, how- ever,only a smallproportionof plots’LAI valueswaswithin the HPDs,i.e.43%forL8,andonly11%forS2.

The LAI retrieval accuracy for the regularizing and informa- tivepriorsaresatisfactory.The datasetweusedinthisstudywas originallyused byKorhonenetal.[20] toempirically predictLAI, amongothervariables, usingaplethora ofvegetationindicesand singleandmultipleband predictors.Theyfound singlebands for bothS2andL8toproducethebestpredictions,withtherededge 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 biasof-0.001.Wefoundthat ourBayesian inversionofPARAS re- sultedinLAIretrievalsofslightlylargerRMSEforL8.Ourbestre- sult, the regularizing prior forL8, 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 Korho- nenetal.[20],whichoriginatefromthefactthatphysically-based models attemptto generalize acrosssensors andviewinggeome- tries, but are not able to account for all influences on the op- tical signal (e.g. atmospheric correction errors) and further need to simplifythe forest radiative transfer process to facilitatecom- putations. The comparison of our results to the empirical study of Korhonen etal. [20] mayindicate the potential of reflectance modelinversionforL8imagesinborealforest.Ontheotherhand, our results for S2 were significantly less accurate than the re- sults of Korhonen et al. [20], both in terms of bias and RMSE.

Besidesthefundamentalmodeldifference betweenourstudyand Korhonenetal.[20],whousedempiricalregressionmodelstore- trieve LAIforthe lidar plots, they also carriedout a band selec- tionwitha penalty termformultiple bandpredictormodels. We didnotaimatselectingindividualbandstotesttheirpotentialfor LAI retrieval accuracy, as this practice seemed to contradict the physical approach, which aims at utilising the entireoptical sig- nal in an attempt to generalize the radiative transfer process in forests.

Inall prior models, Landsat8clearly outperformedSentinel-2 inretrievingLAI, withthebiascomprisingthemajordifferencein accuracy.We attributeS2’s tendencytowardsahigherbiasinLAI retrievals to two characteristics; (i) the differences between the sensors inthe visible bands,and(ii) theincreased spectral reso- lution.Theformer,asdiscussedindetailabove, isduetoS2visi- blereflectancebeinglowerthanL8inourdata.PARASmodelhas shownto have a tendencyto overestimate reflectance ofa given standinborealforests,whichinturnleadtoan overestimationof

leaf area index [15,46]. Hence, a lower reflectance in the visible bandsofS2leadtoanincreasedabsolutebiasintheLAIretrievals.

In the Bayesianinversion approach, all modelparameters are re- trievedsimultaneously. Ourresultsindicatethat theinversionap- proachwasabletopartlycompensate forthelowerreflectancein S2byreducingtheleafalbedointherespectivebandsofconifers, whichwerethedominantspecies.However, asthespectralbands were modeled as partly correlated, the reduction in the visible bandsalsoaffectedtheotherbands,whichmayhaveinhibitedthe compensationforthedifferenceinvisiblebandsofS2.

Modelingthespectralbandsaspartlycorrelatedwasnecessary, asexploratoryanalysisofmodelingthespectralbandsasindepen- dentindicatedanalternativesolutionwhichmistookdeciduousfo- liagewithhigheralbedoforunderstoryvegetation,thusleadingto substantialerrorsintheLAIretrievals.Thesecondpotentialreason fortheincreasedbiasisthehigherspectralresolutionofS2.While ahigherspectralresolutionwouldbeexpectedtoyieldhigherac- curacy of inversion results[37],we observed the opposite to be the case. This issue maybe unique to our dataset andmodeling approach.The increasedspectralresolutionofS2 isexclusivelyin therededgeregion.Thisspectralregionischaracterizedbyasteep riseinleafalbedo,fromlevelsofabout0.1toabout0.8onawave- lengthrangeofabout50to80nm,constitutingthesteepestslope inleafspectraintheopticalregion.Duetothissteepness, anyer- ror,potentiallyoriginatingfromtheinputspectraldata,orthesen- sorspectralresponsefunctions,mayhavehadanexceedinglylarge influenceonthe resultsrelativetoother bandswhicharelocated inrelativelyflatspectralregions.

In addition, we transformed the spectral data to the sensor- specific bandsbefore applyingthePARAS model,whichmayhave hadastronginfluenceduetoits nonlinearity.Thisaprioritrans- formationwasunavoidablefortechnicalreasons,andwasexpected tonotresultinsignificanterrors.However, ashiftintherededge inflictionpointbyafewnanometersmayalreadycausesignificant changesinthespectra,whichwereimpossibletomodelinourap- proach.Wecarriedoutanexploratoryanalysis,whereweadapted thevisiblebandsofS2tothelevels ofL8usingthelinearregres- sionequationsinFig.1aswell aschose asubset ofthebandsof S2thatcorrespondtoL8.Wefoundthatadaptingthevisiblebands of S2 lead to a large bias reduction, and excluding the red edge bandsfurther slightlyreducedthebias,however,atbest thebias ofS2gotclosetothatofL8,butneverless.

Basedonourresults,itappearsthatthesaturationeffectinre- lationtoLAIisadominantfactorthatistobe takenintoaccount when inverting a forest reflectance model. Our recommendation is toregularize LAI wheninverting a forest reflectancemodel, as ourresults show.Theregularizing prior retrievedLAIbest forL8, andonlyslightlyworse inS2.Comparedtousingarealistic prior, which resultedinvery similar accuracyforLAIretrieval,the reg- ularizing prior seems to be favorablefrom a theoretical point of view. When using a regularizing prior, it is not necessary to ob- tain information on a realistic LAIdistribution in the studyarea, butrather it ispossible todefine 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 for- estreflectancemodels.Thisphysicalfoundationofferstheoretically higherpotential forgeneralization becauseitis onlynecessaryto understand when the optical signal saturates, ratherthan having to make assumptions specific to a study site. However, previous work withaBayesian reflectancemodelinversionby Varviaetal.

[46]retrievedLAIusingauniformpriorusingEO1-Hyperionhyper- spectraldata.Hyperiondatahasbeenfoundtoproducelessbiased estimatesofleafbiochemicaltraitsthanLandsat8inaBayesianin- version approachusingPROSPECTmodelatleaflevel[37].Shiklo- manovetal.[37]furtherreportedagenerallyincreasingbiaswith decreasingspectral resolution.Thismayindicatethat whenusing

Fig. 5. Posterior point estimates for conifer share and clumping index.

datawithlowspectralresolution,thechoiceofpriorsbecomesin- creasinglyimportant.

Thereported valuesforLAIretrievals followthe definitionfor the trueLAI, i.e.the unit hemisurface area ofall leavesper unit groundarea.Inmorerecentpublications,theLAIhasbeentermed plant area index, which includes the contribution of other plant elementssuch as woody parts. We acknowledge the importance of accounting forthe contribution of non-leaf plant elements to theforest radiationregime,butusedthetermsLAIandplantarea index interchangeably. The true LAI is the product of the effec- tive LAI [3] and a clumping index. The effective LAI can be in- terpreted astheLAI ofa standwithequivalent angulargapfrac- tions butpurely randomly distributedfoliage, which is necessary for utilizing the modified Beer’s law [35] in the forest radiation regime.Theclumpingindexisthesinglefactordescribingthespa- tial aggregation,or non-randomness of the leaf distribution in a canopy. Hence, a retrieval of true LAI is necessarily a joint pre- dictionofeffectiveLAIandclumpingindex. Thethreeprior mod- elsweusedwerefortheeffectiveLAI,astheclumpingindexwas modeledseparately,andthetrueLAIwasderived fromthosetwo factorsasL=L_e/CI.Therefore,ourresultsforLAIretrievalcanonly be interpreted in combination with the estimates for clumping index.

The posterior modes for both sensors show a similar rela- tionship between the conifer share and clumping index of plots (Fig.5).Therangeoftheclumpingindexstretchesfurtherinboth directions for S2, which also exhibits a slightly wider range of conifersharevalues.Basedonthepriorformulation,therelation- shipbetweenconifershareandclumpingshouldbea linearmix- ture,withaclumpingindexofaround1inpuredeciduousstands, and0.6inpureconiferstands.However,ourresultsshowthatthe clumpingindexdecreases strongerthantheprior formulation,i.e.

ataconifer sharebetween0.7and0.8, thereisadropinclump- ing index down to a posterior mode of slightly below 0.5, with the highestconifer shares in the posteriorat around 0.83. Com- monly, clumping correction is done only for conifer shoots, fol- lowingOker-BlomandSmolander[29],whileeffortstounderstand clumpingofplantelementsathigherhierarchicallevelshavebeen few [42]. While considering the overall modeling uncertainty of thisstudy,theposteriormarginaldistributionoftheclumpingin- dex of results data hints at increased clumping that ranges be- yond shootlevel, i.e.atclumping indicesbelow0.4, whichis the lower endofobserved shoot clumping[29,30].While thesefind- ingsmustbe interpretedwithcautionduetoalackofvalidation, the retrievals of true LAI can be indirectly used to validate the

possibility ofclumping beyondshoot level having a stronger in- fluenceonthe forestradiation regimeasdescribedby thePARAS modelthanpreviouslythought.Furthermore,theseresultsareonly reasonablyobtainable byusing a Bayesian inversionapproach, or anyotherinversionapproachthathasinherentcapabilitiesforun- certaintyquantification.

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 ap- proach, i.e. the spectral properties of leaves and understory be- ing modeled across all stands, did not allow for variation in op- ticalpropertiesbetweenstands.Thismodelingdecisionwasmade duetoboth computationalandmodeling reasons.Modeling opti- cal properties at stand level would have, dueto the large num- ber of plots, required enormous amounts of computer memory thatwere not available. Furthermore,modeling optical properties atstand level would have potentially lead to parameter identifi- cation issues, which we encountered during an exploratory test wherewe modeled all parameters at stand level. The parameter identificationproblemswereparticularlyprominentinthedistinc- tionbetweentreespecies,leadingtoconfusionbetweenconiferous anddeciduousdominatedstands,andsubsequentlytoproblemsin clumpingindexand LAIretrievals. As we were particularlyinter- ested in speciesdistinction in relation to clumping, we chose to modelopticalpropertiesoverthewholestudyarea.Theunderstory opticalpropertieswere modeledaccordingly,leadingtotheinver- sionmodel being unable to distinguish differencesin understory types.This,similarly, wasnecessarytoavoidparameteridentifica- tionproblemsas aconsequence ofthelow spectral resolution of thesensors. However, this modeling choice waspractical for our study,andwepointouttheneedforfurtherinvestigationonaro- busthierarchicalmodelingframeworkforreflectancemodelinver- sion that utilizes synergies betweenindividual stands andwhile efficiently accountingfor and enabling the analysis of variability betweenstands.

Building on the work ofVarvia etal. [45,46], we applied the Bayesianinversionapproachtomultispectraldata.However,adapt- ingthe approachdevelopedforhyperspectral datarequiredusto changetheir prior model.Comparedto Varviaetal. [45],we de- finedpriorsforoptical propertiesthatwerelesscorrelated,which is reasonable given the further distance in wavelength between multispectral and hyperspectral data. We further enhanced the modelto take intoaccount treespecies shares at stand level,to allow a more in-depth analysis of which combinations of forest standparameters likelyconstitute theremotely sensedsignal. On some parameters such as the clumping index, we imposed less constraintsontheprior,i.e.bytruncatingthedistributionatvalues 0.1and1.1. Duetothelower spectral resolution,we modeled the uncertaintyaroundtheobservedbottom-of-atmospherereflectance values with an interval of 20% to account for uncertainty intro- ducedbynoiseanderrorsfromtheimagingprocess,preprocessing oftherawsensordata,theatmosphericcorrection,andthePARAS reflectance model itself. We testeduncertainty parameterizations fortheBoAreflectance between5% and25%in stepsof5%, with 20%showingthebestresultsinterms ofbothnumericalstability andaccuracyofLAIretrievals.

A currentlimitationoftheBayesian approach,particularlythe MCMCmethodweused,theno-U-turnsampler,concernsthesam- plingofmultimodal distributions.There are MCMCmethods that areable to better samplefrom multimodaldistributions, such as inter-chainadaptation[5]orsequentialMonteCarlo[7],however, they were not implemented in the software in a wayusable for ourstudy at the time we carried out the computations. The fo- cus of this study was, in terms of software, to build a Bayesian reflectancemodel inversion approachthat relies onreadily avail-

ablesoftwaretoillustratethepotentialforpracticalapplicationsin remotesensingofvegetation.Wethusrefrained fromimplement- ing custom methods other than implementing the PARAS model intheBayesian framework.The problemofsamplingfrommulti- modaldistributionswithMCMCwasavoidedinthisstudybyusing initiallylooseboundariesofpriordistributions,andempiricallyan- alyzingtheoccuringmultiplemodes,andapplyingdomainknowl- edgeto identifyillogicalsolutions tothereflectance modelinver- sionby settingboundariestothepriors, i.e., truncatingprior dis- tributions atvaluesthatexceed anecologicallymeaningfulrange.

Weexpectthat,withtherapiddevelopmentinBayesiancomputa- tionalmethods andtheir increasing availability toscientific com- munities beyond mathematics and computer science, the prob- lemofmultimodalitywillhopefullybecomelessprominentinthe future.

Finally,therearelimitationsto theresultspresentedherethat originate from the PARAS model itself. The LAI retrieval bias to- wards overestimation as well as the mostly lower leaf and un- derstory reflectance factors indicate that PARAS has a tendency to overestimate reflectance ofa givenstand. There area number ofreasonswhyPARAS mayoverestimate forest reflectance,which comefromthestillpoorunderstandingofsomeparametersinthe model,andthemodelformulationitself.PARASemploysaparam- eter for scattering uncertainty; Q (Eq. (5)), which is based on a semi-physicalrelationshipbetweenscatteringinforwardandback- warddirectionof acomplex medium parameterizedby empirical ray tracing model ina Scots pine stand [28]. While parameteriz- ing the canopy bidirectional scatteringphase function using Q is currently the mostreasonable assumption, thisparameter would greatly benefitfrom further study.Next,it is commonto correct coniferousforest LAIforclumpingatshootlevel,however,clump- ingathigherlevelsmayplayanimportantrole,butissofarpoorly studied [42].Few studies [e.g.24,52] were looking intoclumping basedonterrestriallaserscanning,wehopethattheseeffortswill continueandbecome moreextensive.Clumping, inPARAS model and,moregenerally,inforestradiationregimemodels,playsacru- cial role asit governs therelationship between thetrue LAIand angular gap fractions, which wouldbe the singleparameter, de- spite the leaf angle distribution, that is able to compensate for theinconsistentestimationaccuracybetweenLAIandcanopygaps we found in thisstudy. Finally,though toa lesser extent, PARAS doesnotaccountformultipleinteractionsofphotonsbetweenthe canopy and understory. Multiple interactions have been studied inPARAS andfoundto contributesignificantly tostand albedoin conditions where the understory is covered in snow, but other- wisemultipleinteractionsmaynotplayastrongroleinforestre- flectance[26].

BesidestheseissueswiththePARASmodel,we expect,though to a lesser extent,that theassumptions around using Beer’slaw andtheG-functionforangularleafprojectionarea ratiomayhave contributed to the inconsistencies in the angular gap estimates.

Beer’s law, together withthe G-function, which incorporates the leaf angle distributions, was used to estimate angular gap frac- tions based onLAI. Here again, clumping plays a significant role, asitincreasesordecreasesthegapfractionatagivenLAIandleaf angledistribution.The G-functionwe usedforconifer trees were based onStenberg[39],who derived a version of theG-function for non-flat leaves.While usingthese shapes forconifer needles seems to be more reasonable, the applicationof thisG-function stronglyaffectsthecommonassumptionsofasphericalleafangle distribution,i.e.thatGisnotconstantforanyangle.Furthermore, theconiferG-functionofStenberg[39] has,toourknowledge,not beenappliedinreflectancemodelinversionstudies.Lastly,theleaf angle distributions we used were planophilefor deciduous trees, and sphericalfor conifers. While the planophiledistribution was clearlyrecommendedbyPiseketal.[31]forbroadleaved treesin