Remote Sensing of Environment The SPART model: A soil-plant-atmosphere radiative transfer model for satellite measurements in the solar spectrum

Radiative transfer models (RTMs) of vegetation canopies can be applied for the retrieval of numerical values of vegetation properties from satellite data. For such retrieval, it is necessary ﬁ rst to apply atmospheric correction to translate the top-of-atmosphere (TOA) satellite data into top-of-canopy (TOC) values. This atmospheric cor- rection typically assumes a Lambertian surface re ﬂ ection, which introduces errors if the real surface is non-Lambertian. Furthermore, atmospheric correction requires atmospheric characterization as input, which is not always available. In this study, we present an RTM for soil-plant-atmosphere systems to model TOC and TOA re ﬂ ectance as observed by sensors, and to retrieve vegetation properties directly from TOA re ﬂ ectance skipping the atmosphere correction processes with the inversion mode of the RTM. The model uses three computationally e ﬃ cient RTMs for soil (BSM), vegetation canopies (PROSAIL) and atmosphere (SMAC), respectively. The sub- models are coupled by using the four-stream theory and the adding method. The resulting ‘ Soil-Plant-Atmosphere Radiative Transfer model ’ (SPART) simulates directional TOA spectral observations, with all major e ﬀ ects included, such as sun-observer geometries and non-Lambertian re ﬂ ectance of the land surface. A sensitivity anaylsis of the model shows that neglecting anisotropic re ﬂ ection of the surface in coupling the surface with atmosphere causes considerable errors in TOA re ﬂ ectance. The model was validated by comparing TOC and TOA re ﬂ ectance simulations with those simulated with the atmosphere-included version of the DART RTM model. We show that the di ﬀ erences between DART and SPART are less than 7% for simulating TOC re ﬂ ectance, and are less than 20% (less than 10% at most bands) for simulating TOA re ﬂ ectance. The model performance in retrieving key vegetation and atmospheric properties was evaluted by using a synthetic dataset and a satellite dataset. The inversion mode allows estimating vegetation properties along with atmospheric properties and TOC re ﬂ ectance with reasonable accuracy directly from TOA observations, and remarkable accuracy can be achieved if prior information is used in the model inversion. The model can be used to investigate the sensitivity of surface and atmospheric properties on TOC and TOA re ﬂ ectance and for the simulation of synthetic data of existing and forthcoming satellite missions. More importantly, it facilitates a quantitative use of remote sensing data from satellites directly without the need for atmospheric correction.


Introduction
Earth observation satellites in the optical domain enable monitoring of the status of vegetation canopies. Since the turn of the century, the availability of measurements of reflected solar radiation from satellite platforms has increased spectacularly. We have rapid access to time series observations from various multi-spectral satellite sensors (Slater, 1980;Dorigo et al., 2007), e.g., Landsat since 1972 (Helens and Fires, 1975), Satellite Pour d'Observation de la Terre (SPOT) since the mid-1980s (Chevrel et al., 1981), the Moderate Resolution Imaging influenced by the optical properties of the atmosphere, vegetation and soil or understory background (Verhoef and Bach, 2007). Consequently, to understand TOA signals, not only the radiative transfer in vegetation canopies but also in the atmosphere, as well as soil-canopy and canopyatmosphere interactions, have to be considered. In order to better interpret remote sensing data, a wide range of vegetation and atmospheric RTMs have been developed. A vegetation canopy can be represented simply by a turbid medium (Allen et al., 1970;Suits, 1971;Verhoef, 1984), a two-layer medium (Kuusk, 2001), a multi-layer medium (Wang and Li, 2013;Yang et al., 2017) or explicitly by an 3-D realistic natural scene (Gastellu-Etchegorry et al., 1996;North, 1996). Radiative transfer problems in the canopy can be solved numerically using techniques like Monte Carlo ray tracing (North, 1996;Disney et al., 2000) or analytically using techniques like the four-stream theory (Verhoef, 1985) and discrete ordinates (Knyazikhin et al., 1992). These have resulted a number of vegetation RTMs, such as the Suits model (Suits, 1971), the SAIL model (Verhoef, 1984), the DART model (Gastellu-Etchegorry et al., 1996 and the FLIGHT model (North, 1996). There are also many atmospheric RTMs that link surface signals with TOA signals, e.g., DISORT (Stamnes et al., 1988), SBDART (Ricchiazzi et al., 1998) and LBLRTM (Clough et al., 2005). Among all atmospheric RTMs, MODTRAN (MODerate resolution atmospheric TRANsmission, Berk et al., 2005) and 6S (Second Simulation of a Satellite Signal in the Solar Spectrum, Vermote et al., 1997) have found a wide application in the field of remote sensing due to their reasonable compromise between model simplicity and realism.
The most common way of using vegetation and atmosphere RTMs to retrieve vegetation properties from satellite data is sequential: first atmospheric correction is carried out to obtain a top-of-canopy (TOC) signal by using an atmospheric RTM, followed by the retrieval of vegetation properties by using a vegetation RTM. The problems with this approach are that (i) atmospheric correction requires measurements of atmospheric properties (e.g., aerosol properties and water vapour content) which are unknown in many cases, (ii) most operational atmospheric correction methods assume Lambertian surface reflection and the TOC reflectance obtained after the correction is at best an approximation of the actual directional reflectance of a canopy (Wang et al., 2010). The conventional approaches to acquire atmospheric properties are either based on local or regional atmospheric photometry (e.g., AERONET, AErosol RObotic NETwork, Holben et al., 1998) or satellite data. E.g., the MODIS water vapour products are computed from the near-infrared channels (i.e., band 18 from 915 to 965 nm and band 19 from 915 to 965 nm), and the aerosol products are estimated from the blue channels (Gao and Kaufman, 2003;Vermote et al., 2016). These atmospheric measurements pose a limitation on retrieving surface properties from remote sensing data due to a potential mismatch in both acquisition time and geolocation between the TOA satellite data and the atmospheric properties used for atmosphere correction. Additionally, the Lambertian surface assumption in atmosphere correction has clear drawbacks on retrieving surface properties from remote sensing data as well. It omits the critical information of canopy structure contained in anisotropic reflectance Bach, 2003, 2007;Wang et al., 2010). Although using multi-angular observations (e.g., the Multi-Angle Implementation of Atmospheric Correction algorithm, MAIAC, Lyapustin et al., 2012) allows performing atmospheric correction and retrieving surface bidirectional reflectance distribution function (BRDF) without the Lambertian assumption, the conditions to apply this approach are restrictive: it requires at least four observations of a surface with distinctive sun-observer geometries, and the surface BRDF properties should remain relatively unchanged (Schaaf et al., 2002).
An invertible combined surface and atmospheric model can solve these problems. In such a model, both the surface bidirectional reflectance or BRDF (Nicodemus, 1970;Schaepman-Strub et al., 2006) and the effects of illumination and viewing geometries on the atmospheric paths are considered. The boundary condition of atmosphere radiative transfer is formed by the spectral responses of the components of soil and vegetation. Inverting such a model allows estimating vegetation properties along with the atmospheric properties (Zhang et al., 2009;Verrelst et al., 2015).
The idea of combined vegetation-atmosphere forward modelling is not new. Bach (2003, 2007) coupled SAIL-based reflectance models (i.e., soil-leaf-canopy, SLC) with MODTRAN to simulate TOA observed radiance. Furthermore, MODTRAN output has been integrated into DART in an original atmosphere-earth coupling and mixing model with discretized atmospheric layers and cells, to separate the contribution of atmosphere in ground-based, airborne and spaceborne radiance measurements Yin et al., 2017;Gastellu-Etchegorry et al., 2017;Morrison et al., 2020). These studies have pointed out the potential of these coupled models to retrieve vegetation and atmospheric properties from TOA radiance. However, the use of such coupled models for the retrieval is not yet routine due to the high computational cost of coupled RTMs. Mousivand et al. (2015) and Verhoef et al. (2018) attempted to retrieve vegetation properties from TOA radiance with soil-plant-atmosphere RTMs. However, in these inversions the atmospheric properties that are input to MODTRAN were considered as known. They were used to simulate the optical coefficients (atmospheric transmittance factors), which were then used to translate TOC reflectance to TOA reflectance. In the following retrieval, only the surface properties were obtained by fitting simulated TOA radiance with measured radiance. This procedure is, strictly, not a combined surface-atmosphere retrieval.
The overall objective of this study is to develop an invertible RTM for the combined system of a soil surface, a vegetation layer and an atmosphere layer, that can facilitate the operational retrieval of surface properties from TOA radiance without knowing atmospheric properties. To achieve this objective, we employ the Simplified Method for Atmospheric Correction (SMAC) developed by Rahman and Dedieu (1994). SMAC is a simple atmosphere RTM based on 6S and has comparable accuracy of 6S but operates much faster (Proud et al., 2010). We coupled SMAC with a soil-leaf-canopy RTM using the four-stream radiative transfer theory framework and the adding method (Cooper et al., 1982;Verhoef, 1985). The 'Soil-Plant-Atmosphere Radiative Transfer model (SPART)' is a combination of three existing, computationally efficient models for soil reflectance, canopy and atmosphere radiative transfer, with two innovations compared to the three models individually: (1) the inclusion of non-Lambertian surface reflectance in the atmospheric SMAC model, and (2) a computationally efficient way of coupling the systems by repeatedly adding layers and redefining the background after each addition. In this paper, we present a detailed description of the model and provide preliminary validation and evaluation of the model in simulating TOC and TOA reflectance, and estimating vegetation properties from TOA satellite observations.

Overview of the model structure
A soil-plant-atmosphere system consists of three components: a soil surface, a vegetation layer and an atmosphere layer. The optical properties of each component are simulated with sub-models, which are the BSM (i.e., Brightness-Shape-Moisture) soil reflectance model (Verhoef et al., 2018), the PROSAIL (i.e., PROSPECT+SAIL) canopy RTM (Jacquemoud and Baret, 1990;Verhoef, 1984) and the SMAC (i.e., Simplified Method for Atmospheric Correction) atmosphere RTM (Rahman and Dedieu, 1994). All the sub-models are internally connected and their connections are summarized in Fig. 1. The integrated model simulates both TOC and TOA reflectance of optical satellite sensors in any given viewing direction. The main inputs of the integrated model as listed in Table 1characterizing the soil surface, canopy structure, biophysical properties of leaves within the canopy, atmosphere and sun-observer geometries. The vegetation and atmospheric models are both based on a turbid medium representation, which means that horizontal heterogeneity is excluded.
Only minor modifications of the sub-models BSM and PROSAIL were needed to enable the integrated approach, but SMAC required more substantial revision. The SMAC model assumes that the atmosphere is bounded underneath by a Lambertian surface. In order to enable the inclusion of an anisotropic vegetation canopy layer, we extended the SMAC model to simulate TOA reflectance of a non-Lambertian surface.

Radiative transfer in soil-plant-atmosphere system
Radiative transfer in the soil-plant-atmosphere system is summarized in Fig. 2. Following the classic four-stream theory, the radiative transfer fluxes include the direct solar flux (E s ), downward diffuse flux (E − ), upward diffuse flux (E + ) and upward 'direct' flux in the viewing direction (E o ). For a summary of the notations used in the four-stream theory, please see Table A.1 in Appendix A. We labelled the two turbid media, notably the vegetation and atmosphere layer, as 'layer 1' and 'layer 2', respectively, and the levels 'Top of Soil (TOS)' or 'Bottom of Canopy (BOC)', 'Top of Canopy (TOC)' or 'Bottom of Atmosphere (BOA)', and 'Top of Atmosphere (TOA)' as levels 1, 2 and 3, respectively. Thus 'level 1' separates the soil background from the vegetation layer, and 'level 2' separates the vegetation from the atmosphere, and 'level 3' is the satellite measurement level. For a summary of the notations used in the adding method, please see Table A.2 in Appendix A. Using this numbering for layers and interfaces, the following set of equations describe radiative transfer in the whole system: and the reflectance of the surface at level 1 (R soil or R * (1)) is described by four reflectance factors (i.e., the factors vary with each other due to different types of incident and reflected fluxes) as The interaction of the four streams with each layer is characterized by the layers' scattering matrix, namely R t , R b , T d and T u . They are given by where the subscripts attached to the vectors in the left matrix refer to the direct solar (s) flux, diffuse (d) flux and flux in the observer's (o) direction, and the subscripts attached to the terms in the right matrix denote downward (d), upward (u), top (t) and bottom (b). For example, ρ x1x2 and τ x1x2 (i.e., x1 and x2 are s, d or o) are reflectance and transmittance of the layer for the case of flux x1 to flux x2. R t and R b are the reflectance at top and bottom of the layers, respectively. T d and T u are the downward and upward transmittance, respectively. The surface reflectance of a soil-plant system is composed of (i) reflection of the downward fluxes at the top of the layer R t (1)E d (2) and (ii) the contribution of the background including the multiple reflection between the bottom of the layer and the background T u (1)  multiple reflection. Therefore, the solution of radiative transfer of the system (Eqs. 1a and 1b) for the reflectance of the surface at level 2 (i.e., TOC reflectance factors) is given as where R * (1) is soil reflectance provided by BSM. R t (1), R b (1), T d (1) and T u (1) are the optical properties of the vegetation layer provided by PROSAIL. This solution can also be achieved by introducing Eq. 1a into Eq. 1b and establishing the form E u (2) = R * (2)E d (2), which is explained in detail in Verhoef (1985).
The complete system of soil, vegetation and atmosphere layers can be solved in a similar way, but now we consider the surface of the soilplant system (i.e., at level 2) as the background and add the atmosphere layer on top of this surface. The reflectance of the surface at level 3 (i.e., TOA reflectance factors) is given as where R t (2), R b (2), T d (2) and T u (2) are the optical properties of the atmosphere layer provided by the modified SMAC model.

Directional TOC and TOA reflectance
The reflectances R * (1), R * (2) and R * (3) consist of four components, depending on the direction of the incident and the reflected flux, as in Eq. 3. TOC reflectance by a sensor is the reflectance in the observational direction, which can be expressed as (Nicodemus, 1970): which implies that TOC reflectance is affected by the ratio of diffuse and direct fluxes at the top of the canopy. This ratio is equal to the ratio of atmospheric transmittances τ ss (2) and τ sd (2). Hence, TOC reflectance in the observational direction is expressed as: where τ ss (2) and τ sd (2) are computed by using the modified SMAC model. It is noted that (2) sd sd ss is the fraction of diffuse light at BOA.
The expressions for the reflectance factors of the soil-plant system R so (2) and R do (2) result from Eq. 5: At the top of the atmosphere, the incident radiation consists of the direct solar flux only, and we are only interested in R so , which can be obtained from Eq. 6: In what follows, we briefly introduce the sub-models with a focus on the terms on the right-hand sides of Eqs. 8 and 10. It is worth noting that a simplified version of Eq. 10 can be obtained by assuming, as in the original SMAC model, a Lambertian surface reflectance (i.e., R sd (2) = R dd (2) = R so (2) = R do (2) = R(2)). In that case Eq. 10 is reduced to: Fig. 2. Flux interaction diagram for a soilplant-atmosphere system. (a): light interaction with a atmosphere layer (layer 2) and a vegetation layer (layer 1). (b): light interaction using the four-stream theory. A rectangle stands for incident radiation and a circle for reflected or transmitted radiation. The rectangles and circles are connected by arrows that indicate the direction of flow, and a symbol next to the arrow indicates the corresponding reflectance or transmittance. R and T are the reflectance and transmittance matrices where the superscripts attached to the matrices denote downward (d), upward (u), top (t) and bottom (b) (i.e., R t (1) and R b (1) are the reflectance at top and bottom of the vegetation layer, respectively. T d (2) and T u (2) are the downward and upward transmittance of the atmosphere layer, respectively).
(2) (2) The Brightness-Shape-Moisture (BSM) model simulates the isotropic reflectance for both dry and wet soil, which means R sd (1) = R dd (1) = R so (1) = R do (1) = R soil . It requires soil brightness (B), soil moisture (SM p ), and two spectral-shape related parameters (φ and λ) ( Table 1). This new model has been briefly introduced in Verhoef et al. (2018), but it has not been fully described, especially the modelling of wet soil reflectance.
In BSM, reflectance of a dry soil surface is simulated using three 'basis spectra' (or the so-called global spectral vectors, GSVs) derived from a global soil spectral library by Jiang and Fang (2019). The three vectors are presented in Fig. 3 where a 1 , a 2 and a 3 are fitting coefficients and G 1 , G 2 and G 3 are the global spectral vectors. Although one can fit any given dry soil reflectance spectrum reasonably with the GSVs, these coefficients have a direct relation with the reflectance spectra, but they are not directly linked with soil composition in the same way as PROSPECT reflectance is related to leaf composition. On top of the GSV model, we assume that soil brightness only affects the 'intensity' of soil reflectance and the 'shape' of soil reflectance is determined by other factors, such as roughness, organic matter content and mineralogical composition, and thus we can separate soil brightness effects from spectral shape effects by means of an intensity-shape transformation (Verhoef et al., 2018): where B is soil brightness, which determines the intensity of soil reflectance. φ and λ are related to the other soil properties. This intensityshape transformation is analogous to transforming Cartesian coordinates to spherical coordinates. Hence, the angles φ and λ are analogous to latitude and longitude on the globe.

Soil moisture effects
We use the GSV approach only for simulating dry soil reflectance but not for simulating wet soil reflectance. The soil moisture effects on soil reflectance is described by using a physically-based model, which is based on the water film coating approach (Ångström, 1925). Although a basis spectra for wet soils is also provided in Jiang and Fang (2019), we use the water film coating method in our study for consistency in radiative transfer modelling: The coupling (dry soil plus water) by means of radiative transfer is similar to how we coupled soil, vegetation, and atmosphere.
The water film coating approach treats wet soil as a system of a dry soil layer covered with a thin water layer. In this approach, the reflectance of the combined dry soil plus water film (i.e., wet soil) includes contributions from 1) a first surface interaction (i.e., Fresnel reflection) from the water film coating the soil particles (ρ 12 ), 2) from the background including the multiple reflections between the bottom of the water layer and the background (i.e., (1 (1 − ρ 21 ) + …). The general formula of reflectance of such a system is given as where ρ 12 and ρ 21 are the Fresnel reflectance from one medium to another medium and '1' and '2' refer to air and water, respectively. R bac is the background reflectance, which is slightly higher than a dry soil reflectance spectrum due to the contribution of Fresnel reflection of the water-soil interface (Lekner and Dorf, 1988). exp(−2κ w kΔ) is the twoway transmittance of a water layer, which varies with the optical thickness, where κ w is the water absorption coefficient, Δ is the optical thickness of an elementary water film and k is the number of elementary water films (i.e., k is natural numbers). The water film coating approach has been applied in several soil RTMs, such as the model in Lekner and Dorf (1988), the model of Bach and Mauser (1994) and MARMIT (multilayer radiative transfer model of soil reflectance) of Bablet et al. (2018). A comparison of the soil models mentioned is summarized in Table 2. In both Ångström (1925) and Lekner and Dorf (1988) light absorption in the water layer is ignored (i.e., no exp(−2κ w kΔ) in Eq. 14). In Ångström (1925) and Bablet et al. (2018), dry soil reflectance is used for approximating R bac , while Lekner and Dorf (1988) and Bach and Mauser (1994) consider the effects of the water-soil interface resulting that R bac is slightly higher than a dry soil reflectance spectrum due to the contribution of Fresnel reflection of the water-soil interface (i.e., water-soil Fresnel reflection).
Compared to the above mentioned models, BSM assumes that the water film thicknesses is spatially irregular (i.e., single water film, double water films, etc.), and k follows a Poisson distribution. For water  films of different thicknesses only the transmission loss due to water absorption is modified (i.e., exp(−2κ w kΔ)), since surface reflectance effects are not influenced by the thickness of the film. The fraction of dry soil area, therefore, is P(k = 0). Taken together, we obtain the soil reflectance as: where μ is equal to the expected value of water film thicknesses and is given as a function of soil moisture.
where SM p is the soil moisture volume percentage as the input of BSM and SM c is soil moisture capacity, which is an indicator of a soil's ability to retain water and is 25 (percentage) in the BSM model.

PROSAIL radiative transfer model
PROSAIL is an integration of the leaf RTM PROSPECT and canopy RTM SAIL (Jacquemoud et al., 2009). PROSPECT simulates leaf reflectance and transmittance from 400 nm to 2500 nm. It requires the inputs of the content of leaf pigments, water and dry matter and leaf internal structure (Table 1). PROSPECT is one of the most widely-used leaf RTMs, and it has been developed further since the first publication (Jacquemoud and Baret, 1990). In this study, we used the most recently published version, PROSPECT-D (Féret et al., 2017). This version includes the specific absorption spectra of chlorophylls, carotenoids and anthocyanins.
SAIL up-scales leaf optical properties (i.e., leaf reflectance and transmittance) to canopy optical signals (R t , R b , T d and T u in Eq.1b) by considering the canopy architecture. It also simulates the reflectance factors of the canopy bounded underneath by a soil background using the approach in Eq. 5. In the SPART model, SAILH (Verhoef, 1998) (i.e., SAIL with hotspot effects) is used. It requires the inputs of the leaf inclination distribution and LAI, the ratio of canopy height to leaf size for hotspot effects and sun-observer geometries.
2.6. SMAC atmosphere radiative transfer model 2.6.1. Original SMAC SMAC operates in a broadly similar manner to 6S (Vermote et al., 1997;Proud et al., 2010), requiring the same seven inputs, which are sun-observer geometry parameters (i.e., θ s , θ o and φ so ), aerosol optical thickness at 550 nm (AOT 550 ), ozone content (U O3 ), water water vapour content (U H2O ) and air pressure (P a ). The essence of SMAC compared to 6S is to reduce model complexity and computation time by simplifying the computation of the process variables (e.g., t g , ρ a , S, T(θ s ) and T(θ o )) in the following equation.
where the notations of the original paper of Rahman and Dedieu (1994) are used, and ρ * is the TOA reflectance in the observational direction, t g is two way gaseous transmission, ρ a is atmospheric reflectance, S is the spherical albedo of the atmosphere, and T(θ s ) and T(θ o ) are total atmospheric transmission in the solar and viewing direction, respectively, excluding the effects of gaseous absorption. For a summary of the main notations used in the SMAC mode, please see Table A.3 in Appendix A.
The process variables are computed by using semi-empirical fitting functions in SMAC. The forms of the functions are still based on analytical solutions of the radiative transfer problems, but the coefficients are fitting parameters while they are computed by using physicallybased functions in 6S, which requires more computational resource (see a comparsion between SMAC and 6S, Proud et al., 2010). For example, gaseous transmission at a given band (λ) in SMAC is expressed as: where U (e.g., U O3 and U H2O ) is the vertically integrated absorbing gas amount. a and n are two band-specified coefficients, which are predefined from an analysis of 6S data for the appropriate spectral band. For a given spectral band a and n are constants and are adjusted to outputs of 6S for each of the gases separately. In the calibration of the coefficients, the'US STANDARD' (Krueger and Minzner, 1976) vertical profiles of temperature, pressure and gas concentration were used. Total transmittances (excluding the effects of gaseous absorption) are formulated as a function of AOT 550 and solar zenith angle or viewing zenith angle. T(θ s ) and T(θ o ) at a given band (λ) are given as where a 10 , a 11 , a 12 , a 20 , a 21 and a 22 are band-specified coefficients for a given sensor and a given aerosol type. This transmission takes into account both Rayleigh and aerosol scattering. Aerosol absorption is implicitly taken into account in a 10 and a 20 , depending on the aerosol model (e.g., continental and desert aerosol models). Spherical albedo of the atmosphere (S(λ)) is similarly simplified through the use of predefined constants in an empirical function of AOT 550 and the air pressure relative to standard atmosphere ( P P a a0 ): where a 30 , a 31 , a 32 and a 33 are constants and are adjusted for a given spectral band and for a given aerosol model. Atmospheric reflectance (ρ a (λ)) is the sum of Rayleigh atmospheric reflectance (ρ ar ) and aerosol atmospheric reflectance (ρ aa ), which are given as a function of molecule and aerosols optical properties and sunobserver geometries. The calculation of ρ ar is the same as in 6S: where τ r is molecular optical depth and is a pre-defined constant for a given spectral band of a given sensor derived from 6S outputs. p r (ξ) is the molecular scattering phase function specifying the angular scattering of light by the atmosphere. It is approximated by using the approach in 6S as a function of sun and viewing angles. Aerosol atmospheric reflectance is a function of sun and viewing angles, the asymmetry factor and the single scattering albedo of the atmosphere, of which the latter two are constant for a given aerosol type and a given spectral band. The details of the calculation of the above mentioned process variables are given in Rahman and Dedieu (1994). The coefficients of the semi-empirical fitting functions are calibrated to 6S for each band of given sensors (e.g., LandSAT, MODIS and SPOT) and are stored in a lookup table and read by SMAC at run time as well as by SPART. Therefore, the SPART model is dedicated to satellite sensors and the model is applicable to a wide range of satellite sensors including those listed in Rahman and Dedieu (1994) and many others. However, it is worth noting that the accuracy of SMAC is not sufficient for sun-induced chlorophyll fluorescence study.

Modified SMAC
To consider the anisotropy of TOC reflectance, we first compare the original SMAC formulation of TOA reflectance with a Lambertian surface (Eq. 17) with the formulation using the four-stream theory (Eq. 11). Despite the different notation, they are of identical form. Comparing these two equations, we can express them in similar notation using the following substitutions: (2) Further, to extend the SMAC model for non-Lambertian surface (i.e., from Eq. 11 to Eq. 6), the key is to derive the diffuse (i.e., τ sd (2) and τ do (2)) and direct (i.e., τ ss (2) and τ oo (2)) transmittances separately rather than the sum of them, which is computed in the original SMAC model as t g (θ s , θ o )T(θ s )T(θ o ). We ignore the gaseous transmission for the moment and estimate direct transmittances excluding the effects of gaseous transmission as recommenced in Rahman and Dedieu (1994): (2) exp( /cos( )) (2) exp( /cos( )) where τ is the total optical depth of the atmosphere layer, which is the sum of aerosol and molecule optical depths and is computed in the original SMAC model as: where a 40 , a 41 and a 50 are constants and are adjusted for a given spectral band and for a given aerosol model. Thus, the diffuse transmittances excluding the effects of gaseous transmission are obtained as: Therefore, we can compute TOA reflectance in the observation direction as Note that TOA reflectance can be converted to TOA spectral radiance (L TOA ).
where E a cos θ s /π is extraterrestrial radiance, which is either given as input or estimated according to actual sun-earth distance depending on the day of the year and mean solar exoatmospheric irradiance (Chapter 2.1.3 of Zheng, 2017).

Sensitivity analysis
Simulating TOA and TOC reflectance spectra allows us to test the relative influence of each of the input parameters that control the spectral response. Sun-observer geometries, leaf chlorophyll content and LAI have been identified as the dominating factors of both TOC and TOA reflectance, and aerosol optical thickness has limited impact on TOC reflectance but strongly affects TOA reflectance (Jacquemoud, 1993;Rahman and Dedieu, 1994). Moreover, LAI and leaf chlorophyll content are two essential parameters for vegetation monitoring in ecological and agricultural applications. Therefore, we conducted a sensitivity analysis to examine the effects of these four parameters on TOC and TOA reflectance, namely sun-observer geometies angles (i.e., expressed as BRDF), leaf chlorophyll content C ab , LAI and AOT 550 . In these simulations, the default values of the model input parameters listed in Table 1 were used unless stated otherwise. We took the spectral characteristics of Sentinel-3 OLCI (Ocean and Land Color Instrument) sensor (Nieke et al., 2012) as an example and the corresponding SMAC coefficients for OCLI was used.

BRDF
Both surfaces and the atmosphere contribute to the bi-directionality of TOA reflectance, and SPART is able to simulate both contributions, including the hotspot effects. Examples of BRDF simulation for wavelengths in the visible red at 665 nm and in the near-infrared at 778 nm using the default soil, leaf, canopy and atmosphere parameters are shown in Fig. 4.
Both the polar plot and the plane plot show clear angular and hotspot effects in both TOC and TOA reflectance at the visible and nearinfrared bands (Fig. 4). In the backward-scattering direction (i.e., θ o = 40°, φ so = 0°), the surface (TOC) reflectance at 778 nm is 1.5 times the reflectance in the forward-scattering direction (i.e., θ o = 40°, φ so = 180°). The difference between TOA and TOC reflectance at this band is small, except when the viewing direction approaches the direction of the incident light. In that case TOC reflectance was substantially higher than TOA reflectance due to the hot spot, which was less profound at TOA after interaction with the atmosphere. The atmosphere layer has more significant effects on the visible reflectance than the near-infrared reflectance as indicated by the larger difference between TOC and TOA reflectance at 665 nm than that at 778 nm. TOA reflectance at 665 nm is generally higher than TOC reflectance at this band, and the difference between TOA and TOC reflectance increases with increasing viewing zenith angle. When θ o = 70°, φ so = 180°, P. Yang, et al. Remote Sensing of Environment 247 (2020) 111870 TOA reflectance at 665 nm is 5 times the TOC reflectance due the effect of the atmosphere. The TOA reflectance becomes rather different if we assume that TOC reflectance is isotropic (Fig. 5). Both the magnitude and BRDF of TOA reflectance differ compared to the case in which surface anisotropy is taken into account (Fig. 4). For the NIR (778 nm), the anisotropy caused by the atmosphere is only limited, and TOA reflectance is still nearly isotropic. TOA reflectance at this wavelength ranges from 0.43 to 0.50, whereas it ranges from 0.34 to 0.50 when surface anisotropy is taken into account. The isotropic assumption for TOC reflectance leads to a relative error in the TOA reflectance as large as 15%. For the red (665 nm), the atmosphere introduces considerable anisotropy, such that strong BRDF effects are present even for a Lambertian surface (Fig. 5). At this wavelength, the atmosphere introduces a strongly angular dependent effect on TOA reflectance. At both wavelengths (665 nm and 778 nm), the hot spot effects disappear when surface anisotropy is removed.
3.2. Sensitivity to chlorophyll content, LAI and aerosol optical thickness Fig. 6 shows that leaf chlorophyll content has almost no effects on TOC and TOA reflectance in the spectral region from 400 to 500 nm and from 800 to 1100 nm, but strongly affects the reflectance from 500 to 800 nm. In this spectral region, both TOC and TOA reflectance decrease with increasing chlorophyll content.
Canopy LAI affects the reflectance in the whole spectral region from 400 to 1100 nm (Fig. 7). In contrast to leaf chlorophyll content, LAI strongly affects near-infrared TOC and TOA reflectance. A change in LAI from 0.5 to 6 causes an increase in TOC and TOA reflectance at 778 nm of about 0.17, and a decrease in TOC and TOA in the visible region. Figs. 6 and 7 confirm that the difference between TOA and TOC reflectance is minor in the near-infrared region, except in the absorption bands of atmospheric gases, e.g., absorption by water centred at 970 nm and by oxygen centred at 761 nm.
The aerosol optical thickness at 550 nm (AOT 550 ) has a minor effect on TOC reflectance, but a large effect on TOA reflectance (Fig. 8). Note that for a Lambertian surface, AOT 550 does not affect TOC reflectance at all, since the minor effect on AOT 550 on TOC reflectance is solely a directional effect via the influence of AOT 550 on the direct to diffuse radiation ratio at BOA (Eq. 10).
TOA reflectance decreases with increasing AOT 550 in the near-infrared region, but increases in the visible region. The effects of increasing AOT 550 on TOA reflectance are similar to the effects of decreasing LAI (Fig. 7), but again, LAI has a larger impact in the nearinfrared region than in the visible region while AOT 550 has equal if not larger impact in the visible region.

Model validation and evaluation
Validation of the SPART model was performed by conducting a model intercomparison experiment. In this experiment, we compared TOC and TOA reflectance simulation results from SPART with those from the DART model for a wide range of synthetic scenarios. In the spectral domain from visible to thermal infrared, DART is one of the most detailed RTMs, which combines unlimited discrete ordinates (streams) with exact kernels (Yin et al., 2013) to produce precise forward modelling of the canopy (Widlowski et al., 2007). Moreover, the atmosphere-surface coupling of DART has been validated against MODTRAN in basic atmospheric configurations , and against actual measurements in an urban heterogeneous scenario (Morrison et al., 2020). Hereafter, the version that includes the atmospheric RTM ) is referred to as DART.
Evaluation of the model performance in estimating vegetation and atmospheric properties was conducted by both using synthetic and satellite TOA measurements. The synthetic TOA radiance data was simulated by using DART and the same vegetation-atmosphere scenarios were used as those in the validation experiment. Because the scenarios were pre-defined and the 'true' values of the model parameters were available, we evaluated the model by comparing the estimated values of the parameters to the assigned 'true' values in the generation of the synthetic dataset. We also applied the model to retrieve vegetation properties from MODIS TOA radiance observations for several study sites, and the estimated vegetation properties were compared with the existing MODIS land surface products.

Synthetic scenarios and DART simulation
A synthetic dataset comprising vegetation and atmospheric properties, sun-observer geometries, and the corresponding reflectance and radiance at both top of canopy and top of atmosphere, was generated. We defined a number of synthetic vegetation-atmosphere scenarios that vary from each other by either leaf chlorophyll content, canopy LAI or aerosol optical thickness. In total, there were 144 scenarios covering all the possible combinations in Table 3. The remainder of leaf, canopy and atmosphere properties were kept as the values in Table 1. Once the scenes were characterized, we employed the DART model to simulate both reflectance and radiance at top of canopy and top of atmosphere using the default sun-observer geometry (i.e., θ s = 40°and θ o = 0°).
Besides varying the vegetation and atmospheric properties, the BRDFs at top of canopy and top of atmosphere of the default vegetationatmosphere scene (Table 1) were simulated with DART as well. This aims at validating the SPART model in simulating BRDFs. For DART simulation, 2100 spectral bands (from 400 to 2500nm) and 60 discrete directions (Yin et al., 2013) were used to model the propagation and scattering of each emitted ray within atmosphere and the vegetation canopy. 30 vertical cells of vegetation canopy were constructed to ensure the leaf area index represented by each cell is less than 0.2. The atmosphere was vertically discredited into layers of 500 m/2000 m depth below/above 4000 m altitude to precisely model the distribution of gases and aerosols. The simulations used the default USSTD76 gas model and RURAL aerosol model. The optical properties of gasses and aerosols were weighted by the solar constant at 1cm −1 spectral interval. The optical properties of vegetation were weighted by the solar irradiance at TOC level. The 2-way atmospheric transmittance was considered and weighted by the coupled radiance at the TOC layer. Broad bands were simulated by integrating the results with spectral response functions. Because DART does not include the BSM soil reflectance model, the background reflectance in the DART simulation was taken from the soil reflectance simulated by using the BSM model  with the default inputs (Table 1), ensuring the same optical properties of the backgrounds. The leaf RTM in DART is PROSPECT-D, which is the same as in SPART.
In addition to the purpose of validating the SPART model in simulating TOA and TOC reflectance, this synthetic dataset was used to evaluate the model performance in estimating vegetation and atmospheric properties as well. In this experiment, we added artificial measurement errors to the DART simulated TOA radiance to obtain realistic noise-contaminated TOA radiance observation of OLCI. A generic sensor noise model that considers the noise variance is as a linear function of the detected radiance was used (Verhoef et al., 2018): where L TOA is observed TOA radiance, and a and b are constants that are independent of any surface or atmospheric properties and depend only on technical sensor design parameters. These constants of OLCI were approximated as a = 1 × 10 −4 and b = 2.5 × 10 −4 noting the units of L TOA and a are mWm −2 sr −1 nm −1 and the unit of b is the square of mWm −2 sr −1 nm −1 following the recommendation of Verhoef et al. (2018).

Satellite measurements and products
TOA radiance observations from MODIS (level-1 products, MOD021KM) of several study sites were used to test the potential of the SPART model for estimating surface and atmospheric properties. The study sites were chosen within the FLUXNET site list (Baldocchi et al., 2001), which represent three plant function types (PFTs), namely evergreen forest (EF), mixed forest (MF) and savanna (SAV), respectively (Table 4). The FR-Pue site is located on a flat plateau with elevation of 270 m located 35 km NW of Montpellier (France). The vegetation around this site is largely dominated by a dense overstorey of evergreen oak (i.e.,Quercus ilex L.), and the mean tree height was about 5 m and overstorey LAI was 2.8 ± 0.4 according to the field measurements by Allard et al. (2008). The ES-LMa site is located in a Mediterranean tree-grass ecosystem in southwestern Spain. This is a managed savannah consisting of sparse trees and a herbaceous layer. Trees (mainly Quercus ilex L.) present a fractional cover~20% and tree distance~18.8 ± 5 m (Pacheco-Labrador et al., 2016). The herbaceous is dominated by the three annual plants, namely grasses, forbs and legumes, which are green and active from October to end of May (Luo et al., 2018;Martini et al., 2019). The CH-Lae site is located at the south slope of the Lägeren with elevation of 682 m about 20 km NW of Zurich (Switzerland). The vegetation around this site is a diverse mixed deciduous and coniferous forest, dominated by Fagus sylvatica L., Picea abies (L.) Karst., Fraxinus excelsior L., and Acer pseudoplatanus L (Heim et al., 2009).
The MODIS data used in this study covers the period from January-2016 to December-2018. MOD021KM (version 6) TOA radiance values of the 20 reflected solar bands (i.e., band 1-19 and band 26) were extracted. We took TOA radiance products centred spatially at the location of the selected sites and excluded the data which are identified as low quality using quality flags included in MODIS level-1 products. Note that the air pressure was estimated from sea level air pressure and the elevation of the site, which are available as auxiliary data in MODIS TOA radiance products, and sun-observer geometries are also available in the MODIS level-1 products.
The retrieved LAI and AOT 550 were benchmarked with MODIS LAI products MCD15A3H version 6 (Myneni et al., 2002) and AOT products MOD04L2 (Remer et al., 2005), respectively. The MCD15A3H LAI product is a 4-day composite data with 500 m pixel size. The operational algorithm of this MODIS product is based on atmosphere corrected TOC reflectance. A look-up-table (LUT) method is used to achieve inversion of the radiative transfer problem, and a back-up method based on empirical relations between the normalized difference vegetation index (NDVI) and LAI, together with a biome classification map, is used when the LUT method fails to localize a solution   (Knyazikhin, 1999;Myneni et al., 2002). The daily MOD04L2 was produced with a spatial resolution of 10 × 10 km at nadir and contains AOT values at three different wavelengths (i.e., 470 nm, 550 nm and 660 nm). The MODIS radiance data (MOD021KM) and MODIS LAI products (MCD15A3H) and MODIS AOT products (MOD04L2) are freely available on NASA (National Aeronautics and Space Administration) data repositories.

Retrieval algorithms
We used the numerical optimization method to retrieve surface and atmospheric parameters from TOA radiance spectra. We examined three different configurations of the numerical optimization method to conduct the retrievals from the DART simulated TOA radiance data, while in the application to the MODIS TOA radiance data, we only applied the most practical and realistic configuration.
The numerical optimization method aims at minimizing a cost function, which quantifies the differences between observed (or DART simulated) (L obs ) and modelled (L mod ) radiance by successive changes of the input parameters that is expressed as A cost function can also include the difference between a priori and estimated values of the model parameters. Prior information about the surface and atmospheric properties can reduce the ill-posedness of the numerical optimization and increase the accuracy of the retrieval, and it is usually necessary to include prior information (Laurent et al., 2014). In order to evaluate the effects of prior on the accuracy of estimating C ab , LAI and AOT 550 , we conducted the retrievals for three configurations.
In the first configuration, the values for all the other parameters besides C ab , LAI and AOT 550 were also unknown. In contrast, the second configuration all model parameters were known by their exact values, specified as prior, except for C ab , LAI and AOT 550 . These two configurations represent two extreme cases: a fully unconstrained and an almost perfectly constrained and retrieval. Reality will be somewhere between these extremes. The third configuration represents such case: The values for all parameters were unknown (as in the first configuration), but prior information dictates that the values for each parameter followed a uniform distribution over the interval defined by the lower (P LB ) and upper (P UB ) boundaries in Table. 1. The standard deviation (σ P ) of a uniformly distributed variable is ≈ 1/ 12 0.3 of the intervals of the parameter. This prior has been commonly used for retrievals from TOA radiance (Verhoef et al., 2018) and TOC reflectance (Van der Tol et al., 2016;Celesti et al., 2018). In order to use this prior, we employed a Bayesian framework by adding elements to the cost function that describe the difference between the current set of variables and the prior expected, normalized by the standard deviations expressing the associated uncertainties.
where σ L is the uncertainties of TOA radiance (Eq. 28). P is the posterior and P 0 is the prior array of parameter values (also the default values).
Note that compared to TOA reflectance, the main advantage of using TOA radiance for retrieving surface and atmospheric properties is that uncertainties of TOA observations can be included in the cost function.
To retrieve surface and atmospheric properties from MODIS TOA radiance measurements, we used the third configuration (Eq. 30) as recommended in Verhoef et al. (2018). To perform the numerical optimization, we used the function 'lsqnonlin' of the optimization toolbox of Matlab R2017a, selecting a Trust Region algorithm for updating parameter values within the ranges shown in Table 1. It operates in several steps: 1) taking an initial guess of each parameters (e.g., the default values in Table 1) and running the SPART model to simulate TOA radiance, 2) computing the value of the cost functions (Eq. 29 or 30), 3) computing the Jacobian matrix of the model, which indicating the amount and direction (i.e., increase or decrease) of a change in the output of the cost function caused by a small change in each model parameter, 4) increasing or decreasing the values of the model parameters according the Jacobian matrix to minimize the cost functions and 5) repeating the previous steps until the improvement of the cost function was less than 10 −3 .

Reflectance simulation results from DART and SPART
As explained in Section 4.1.1, DART and SPART use the same leaf RTM (PROSPECT-D), and the same soil reflectance and very similar canopy characterizations. The TOC reflectance simulated with the two models are very similar (Fig. 9). The absolute difference in the TOC reflectance is less than 0.01 and the relative error is less than 7% across all the 21 bands for all the 144 scenarios, and the absolute difference exhibits a wavelength dependence, increasing towards the near-infrared region.
Compared to the difference in the simulated TOC reflectance, the TOA reflectance simulated by DART and SPART show larger differences (Fig. 10), in particular at the blue bands (around 400 nm) and the oxygen absorption feature (around 760 nm). The divergence of SPART from DART in these bands is evident in Fig. 10b, which shows a scatter plot TOA reflectance of the two models, with the blue bands highlighted in blue color and the the band close to the oxygen absorption feature in red color. Nevertheless, the absolute error is less than 0.01 across all nineteen other bands. It is noteworthy that the large relative error at the red edge (around~690 nm, Fig. 10d) is mainly due to the lowest absolute reflectance (less than 0.1 for all the scenarios) at this band.
The effects of viewing angles on the TOC reflectance in the principal panel are identically modelled by SPART and DART as shown on Fig. 11a and b (for a statistical comparison, see Fig. S1 in the Supplementary materials). The differences between the two sets of TOC reflectance simulation are less than 0.01 at each band for various viewing angles, including in the hot spot direction (θ o = 40°and φ so = 0°). The difference in TOA reflectance from DART and SPART is less than 0.02 (Fig. 11b, d). Generally, the large differences of TOA reflectance (> 0.01) occur when the viewing angle is greater than 40 degrees. The difference is relatively larger in the visible region than that in the nearinfrared region.
5.2. Results of model inversion using synthetic TOA radiance data 5.2.1. Estimated parameters Fig. 12 (a1, b1 and c1) shows that the model can estimate C ab , LAI and AOT 550 well if the other model parameters are known. The relative errors are always smaller than 10% for AOT 550 , 15% for LAI and 25% for AOT 550 . Although the impact of LAI and AOT 550 affects TOA reflectance in a similar manner (decreasing LAI or increasing AOT 550 results into decreases of near-infrared TOA reflectance and increases in visible TOA reflectance) (Figs. 7 and 8), their effects on TOA reflectance are still distinct enough to allow for their accurate estimation from TOA radiance. However, if the other parameters are free as well, the accuracy of estimating C ab , LAI and AOT 550 becomes significantly lower, in particular for C ab and LAI (Fig. 12 a2, b2 and c2). The relative errors in the estimated AOT 550 are still less than 40%, but the relative errors in the estimated LAI and C ab are as high as 65% and 70%, respectively.  Table 3: absolute (a) and relative differences (c) changing with wavelengths; a scatter plot between them (b); and the spectra of one representative scenario where LAI = 3, C ab = 50 μg cm −2 and AOT 550 = 0.9.  Table 3: absolute (a) and relative difference (c) changing with wavelengths; a scatter plot (b) between them and the spectra of one representative scenario where LAI = 3, C ab = 50 μg cm −2 and AOT 550 = 0.9. Note: in the panel b, the reflectance at the blue bands and at the oxygen absorption feature are marked in blue and red, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) P. Yang, et al. Remote Sensing of Environment 247 (2020) 111870 Assuming uniform a priori distributions for the variables over the interval, we estimate LAI, C ab and AOT 550 reasonably well (Fig. 12 a3, b3 and c3). The improvement in C ab compared with the fully unconstrained retrieval is substantial (Fig. 12 a2, b2 and c2). The estimation of AOT 550 is reasonably accurate in all three conditions with relative errors less than 40%. The numerical optimization approach generally underestimates the LAI when LAI ≥4. The use a prior information significantly reduces this underestimation compared to that of the fully unconstrained retrieval (Fig. 12 a2), while the underestimate is the least in the fully constrained retrieval (Fig. 12).

Estimated TOC reflectance
TOC reflectance estimated by using the model inversion approach with the SPART model is further compared to that obtained by using traditional atmospheric correction with the SMAC model, and both approaches are compared to the 'true' TOC reflectance of the 144 synthetic scenarios (Fig. 13). The model parameters are assumed to uniformly distributed in their interval for the model inversion approach in Fig. 13, and the results from other two conditions (i.e., all the model parameters are free and only C ab , LAI and AOT 550 are unknown and the other model parameters are known as prior) can be found in the Supplementary materials (Fig. S2).
Although the atmospheric characteristics (i.e., AOT 550 , U O3 and U H2O ) are well defined in the atmospheric correction, there is still a considerable error in the estimated TOC reflectance due to ignoring the BRDF effects of the TOC reflectance during the correction, in particular in the NIR region. The atmospheric correction approach generally overestimates the directional TOC NIR reflectance of the 144 scenarios and the error is as high as 0.03. There are higher relative errors in the visible bands than the NIR band, but mainly due to the low reflectance in the visible bands. In contrast, the model inversion approach reproduces the TOC reflectance better than the atmospheric correction approach regardless of whether prior information is applied or not (see the results for different inversion approaches in Fig. S2 in the Supplementary materials). Although there are considerable errors in the estimated vegetation parameters (Fig. 12 a1, a2, a3, c1, c2 and c3), the model inversion approach reproduces the TOC reflectance well with an error less than 0.01. If only C ab , LAI and AOT 550 are unknown and the other model parameters are known as prior (model inversion 2), the TOC reflectance estimated by using the model inversion approach is almost identical to the 'true' TOC reflectance in the synthetic dataset with errors less than 0.001. If all the model parameters are free (model inversion 1) or the model parameters are assumed to uniformly distributed in their interval (model inversion 3), the performance is still satisfactory and better than using the atmospheric correction approach.

Results of model inversion using satellite measurements
The estimated LAI from MODIS TOA radiance using the SPART model presents similar seasonality with the MODIS MCD15A3H LAI products (Fig. 14). The MF site has the most pronounced seasonality with low LAI in the winter (from December to March) and high LAI in the summer (from May to September). In contrast, in the SAV site, both MODIS MCD15A3H LAI and retrieved LAI from TOA radiance show low values in the summer. There is much less seasonal variation of LAI in the EF site, where the mean LAI is 2.4 and 2.7 for MODIS-SPART and the MODIS MCD15A3H product, respectively. Additionally, in both MODIS MCD15A3H LAI and SPART retrieved LAI, there are unrealistic fluctuations. The two LAI products are strongly correlated among these three sites with overall R 2 of 0.66 (Fig. 15). In the SVA and EF sites, they are moderately correlated with R 2 of around 0.43, and there is a small seasonal variation of LAI in these two sites (Figs. 14 and 15), while the R 2 is as high as 0.82 in the MF site.
Together with surface properties, atmospheric properties were estimated from TOA radiance. The estimated AOT 550 from MODIS TOA radiance using the SPART model has comparable magnitude (from 0 to 0.5) with the MODIS AOT products (Fig. 16). In both the SVA and EF sites, AOT has a pronounced seasonality with low AOT from September to March and high AOT from March to September. The seasonal variation in these two sites are similar. The two AOT estimates are positively correlated among these three sites with overall R 2 of 0.54 (Fig. 17). In the SVA site, they are strongly correlated with R 2 of around 0.70.  (Table 1) are presented. P. Yang, et al. Remote Sensing of Environment 247 (2020) 111870 6. Discussion

The four-stream theory and adding method for coupling RTMs
The four-stream theory and the adding method provide a practical framework for coupling the soil reflectance model BSM with the canopy RTM PROSAIL and further with the modified atmosphere RTM. A soilvegetation-atmosphere ensemble can be treated as a multi-layer system and the radiative transfer problem in such a system can be solved by using the adding method and four-stream theory. In the SPART model, the soil is considered as a Lambertian surface, and canopy and atmosphere are considered as turbid media for simplicity, but the framework for coupling RTMs can also be easily applied to non-Lambertian soil and to multi-layer canopies or atmosphere (Yang et al., 2017). Furthermore, this framework is applicable to land covers other than vegetation, such as soil and water body (Salama and Verhoef, 2015). An equivalent model for a soil-water-atmosphere system can be achieved by removing the vegetation layer or by replacing the vegetation layer with a water layer.
In order to employ the four-stream theory to characterize radiative transfer in the atmosphere, the SMAC model is extended. In the modified SMAC, diffuse and direct transmittances rather than their sum (total transmittances) are computed separately. This extension of SMAC allows including the bidirectional effects of TOC reflectance in simulating TOA reflectance, which provides additional vegetation information, especially on canopy architecture (e.g., LAI and leaf angle). Neglecting the BRDF effects of vegetation canopies in coupling the surface with atmosphere causes a considerable error in TOA reflectance (Figs. 4 and 5). Moreover, neglecting BRDF effects and assuming a Lambertian surface in atmospheric correction cause errors in TOC directional reflectance, even though we could know atmospheric properties at the location during satellite data acquisition (Fig. 13).

Model applications
Forward simulation of the coupled surface and atmosphere RTMs can be used to (1) evaluate the effects of sun-target-observer geometries, surface and atmosphere properties on TOA and TOC reflectance (as in Figs. 4 and 5), and (2) to analyse the sensitivity of data of satellite missions with realistic synthetic data (as in Figs. 6, 7 and 8). The SPART model offers a very high level of accuracy of TOC reflectance (a relative error of less than 7%, Fig. 10) and a reasonably accurate simulation of TOA reflectance (a relative error of 0 − 20%, Fig. 10), when compared to the DART model. The use of different atmosphere RTMs (i.e., SMAC in SPART and a MODTRAN-based RTM in DART) is the main cause of the difference in TOA reflectance. The results agree with the findings of a comparison between SMAC and 6S, which shows a up to 20% difference in atmospherically corrected TOC reflectance in many sets of atmospheric conditions due to some simplifications in the SMAC model (Proud et al., 2010). Additionally, both the present study and Proud et al. (2010) find that the simulations involving medium-to-high viewing zenith angles (greater than 40°) contain a higher relative error than the simulations at lower viewing zenith angles. Despite fact that SMAC uses semi-empirical fitting of 6S output with pre-calculated  Table. 1 are known except LAI, C ab and AOT 550 . Middle panels: all parameters in Table 1 are unknown. Lower panels: all parameters in Table 1 are unknown but assuming uniform a priori distributions for the variables over the interval. Note that in all three conditions the sun-observer geometries are known.
coefficients for specific sensors, the performance is still satisfactory as demonstrated by a comparison to 6S (Proud et al., 2010) and the results from Section 5.1. The main advantage is the three orders of magnitude lower computational demand than 6S, as shown in Rahman and Dedieu (1994).
The SPART model has obvious advantages in the inversion mode Fig. 13. Estimation of TOC reflectance from synthetic Sentinel-3 OLCI TOA reflectance by using the atmospheric correction approach and by using the model inversion approach.   Yang, et al. Remote Sensing of Environment 247 (2020) 111870 due to its simplicity. Three sub-models and the SPART model itself require low computational resources. BSM and PROSAIL are analytical models and their simplicity facilities the retrieval of vegetation properties from ground and airborne TOC reflectance measurements (Celesti et al., 2018;Yang et al., 2019). Even without further optimization of the numerical inversion algorithm, it requires only a few seconds to retrieve one set of parameters reasonably accurate from a reflectance/radiance measurement on a regular PC. Besides the high computational efficiency, the vegetation-atmosphere coupled model allows retrievals from TOA observations directly without atmospheric correction. The combined use of TOA observations and SPART provides similar seasonality and comparable magnitudes of the MODIS LAI products (Figs. 14 and 15), and AOT products (Figs. 16 and 17). From this perspective, the SPART model can be used for (1) retrieving vegetation and atmosphere properties from existing satellite data, (2) optimizing the specifications of sensors and data acquisition scenarios and (3) experimenting with various retrieval algorithms and data assimilation procedures.
SPART takes sensor characteristics and observational geometries quantitatively into account, by respectively spectral convolution and radiative transfer modelling. This makes observations from different sensors and observational conditions comparable. Consequently, timeseries observations from multiple sensors and orbits can be combined relatively easily. This may help reduce the problem of ill-posedness in the retrieval process. The accuracy of the retrieval not only depends on the number of free parameters (Fig. 12), but also on the number of observations of a target, where a larger number of observations of the same target reduces the degrees of freedom.
Additionally, the coupled model facilitates the use of prior information, which is expected to improve the accuracy in estimating the parameters of interest. This is particularly useful in retrieval from timeseries observations. In time series of consecutive remote sensing measurements, prior information may be available from retrievals in previous time steps. Earlier retrieved parameter values may help constrain the inversion problem at subsequent time steps. Here one can make use of the more gradual change in land surface properties compared to atmospheric properties, by assigning resistances to change between observations of the same location. This may help to improve retrieving vegetation properties from time-series satellite measurements and reduce unrealistic fluctuations (Fig. 14), which exist in most of the satellite products. Prior information from time series retrieval may complement more general prior information based on land cover classification. For example, forest tends to have high LAI and croplands have clear seasonal cycle of both leaf chlorophyll content and LAI.
The simplicity of the SPART model has many advantages, but also several limitations. Both vegetation and atmosphere are considered as turbid media, and the model does not consider the variation in other atmospheric constituents than ozone and water. In this study we performed two limited model inversion exercises, one using synthetic data and one using MODIS TOA data. The limitation of using synthetic data is that the dataset only covers the scenarios of homogenous canopies and cloud-free conditions. Investigation of the accuracy of the model for more complex canopies and atmospheric conditions is required. Nevertheless, it is noteworthy that 1D and 3D models tend to converge at coarser spatial resolutions (typically beyond 100 m) and the simple PROSAIL model appears suitable for satellite applications due to its simplicity and reasonable accuracy (Widlowski et al., 2005). The retrieval from MODIS data that we carried out was limited in the sense that independent validation data were lacking. The MODIS LAI products can only serve as a benchmark, but not as the true values since they are derived from the same MODIS observations by model inversion. For example, the MODIS LAI is derived from remote sensing signals as well (but with different approaches), which may suffer from the same 'saturation issue' when the LAI is high (e.g. > 4). The retrievals in this study were conducted by using the algorithm that only uses uninformed prior information to constrain the retrievals. Further work in this direction should include the development of an algorithm for the retrieval of specific parameters for specific land cover.

Conclusions
This study presents a simple radiative transfer model to simulate the effects of soil background, leaf properties, canopy structure, atmosphere properties and sun-observer geometries on TOA and TOC reflectance observed by sensors. The benefits of the model presented are simplicity and improved treatment of the interaction between atmosphere and canopy. The resulting SPART model is computationally efficient and easy to invert.
Compared to other coupled models published before, such as the combined PROSAIL and MODTRAN model and the DART model with an atmosphere layer, the computational effort is orders of magnitude less while maintaining the essence. This computational efficiency can be attributed to not only the three sub-models, but also the framework of using the four-stream theory and the adding method for coupling the sub-models. The SPART model may be used in operational inversion schemes to retrieve the surface and atmospheric properties from timeseries data from multiple sensors without atmospheric correction.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. direct transmittance in the direction of the sun τ sd diffuse transmittance for specular incidence τ dd diffuse transmittance for diffuse incidence τ do directional transmittance for diffuse incidence τ oo direct transmittance in the direction of observation ρ so bidirectional reflectance of a layer ρ sd diffuse reflectance for specular incidence of a layer ρ dd diffuse reflectance for diffuse incidence of a layer ρ do directional reflectance for diffuse incidence of a layer R so bidirectional reflectance of a surface R sd diffuse reflectance for specular incidence of a surface R dd diffuse reflectance for diffuse incidence of a surface R do directional reflectance for diffuse incidence of a surface   Terminology used in the SMAC model (after Rahman and Dedieu (1994)).