A Refined Four-Stream Radiative Transfer Model for Row-Planted Crops

In modeling the canopy reflectance of row-planted crops, neglecting horizontal radiative transfer may lead to an inaccurate representation of vegetation energy balance and further cause uncertainty in the simulation of canopy reflectance at larger viewing zenith angles. To reduce this systematic deviation, here we refined the four-stream radiative transfer equations by considering horizontal radiation through the lateral “walls”, considered the radiative transfer between rows, then proposed a modified four-stream (MFS) radiative transfer model using single and multiple scattering. We validated the MFS model using both computer simulations and in situ measurements, and found that the MFS model can be used to simulate crop canopy reflectance at different growth stages with an accuracy comparable to the computer simulations (RMSE < 0.002 in the red band, RMSE < 0.019 in NIR band). Moreover, the MFS model can be successfully used to simulate the reflectance of continuous (RMSE = 0.012) and row crop canopies (RMSE < 0.023), and therefore addressed the large viewing zenith angle problems in the previous row model based on four-stream radiative transfer equations. Our results demonstrate that horizontal radiation is an important factor that needs to be considered in modeling the canopy reflectance of row-planted crops. Hence, the refined four-stream radiative transfer model is applicable to the real world.


Introduction
Canopy reflectance modeling for row-planted crops is a hot topic in quantitative remote sensing of vegetation, which is promising for developing physics-based methods to estimate canopy biophysical parameters [1][2][3][4]. In the inversion of remote sensing, the efficiency of the inversion is closely related to the calculation time and accuracy of the model [5]. In general, the physical modeling methods can be separated into three categories, namely geometric optical, computer simulation, and radiative transfer methods [5]. In these modeling techniques, compared with geometric optics that only consider surface reflectance [6,7], radiative transfer considers the leaf inclination angle distribution function (LADF) from a statistical perspective inside the canopy [8,9]. Then, based on LADF, the volume scattering inside the canopy can be accurately described [10], thereby improving the accuracy of the reflectance. Simultaneously, compared with computer simulation, radiative transfer considers the statistical LADF inside the canopy, then the volume scattering can be described inside the canopy. Although radiative transfer based on volume scattering from a statistical perspective has a lower 2. Methodology

Four-Stream Radiative Transfer Equations for Continuous Crops
The four-stream radiative transfer theory was developed based on the K-M theory (it is the equations proposed by Kubelka, P. and Munk, F.; hence, it is named as K-M equations) [16,50]. The theory assumes that the canopy consists of an isolated homogeneous scattering layer (i.e., plane-parallel medium), which is composed of small and flat leaves (Figure 1c). The four-stream radiative transfer equations are written as [13]: Remote Sens. 2020, 12, 1290 4 of 28 Here, L is the leaf area index (LAI) and   l f  is the leaf inclination distribution function (LADF), which uses an ellipse distribution [51] (This is an equation for azimuthally independent analysis [19]).
Since s E , E  , E  , and o E were assumed to be vertical in the four-stream radiative transfer equations (Z-axis in Figure 1(d)), Equation (1) describes the extinction of specular flux in the vertical direction, Equations (2) and (3) describes the vertical transfer process of the upward and downward diffuse fluxes, and Equation (4) is an approximation of the radiative transfer equation, referring to the entire vertical radiative transfer process from the viewing [52]. Combined with the assumption that the canopy was infinitely extended (Figure 1(c)) for four-stream radiative transfer equations, radiative transfer of vertical flux density achieved good results in continuous crops. However, for the row canopy, there is a radiative transfer parallel to the horizontal ground surface and passing through the lateral "walls" (Figure 1(a)). Therefore, the lack of horizontal radiative transfer may cause an energy imbalance when crops have a row structure, which affects the accuracy of the model. Here, A 1 is the row width, A 2 is the distance between rows, and h is the height of the canopy.
Here, k and K are extinction coefficients; s , s, σ, w, v, and v are scattering coefficients; and a is the attenuation coefficient. These parameters are optical coefficients in Equations (1)-(4) [19]. For unknowns in Equations (1)-(4), E s , E − , E + , and E o are the downward specular irradiance, downward hemispherical diffuse flux density, upward hemispherical diffuse flux density, and flux-equivalent radiance in the viewing direction, respectively. Finally, L is the differential leaf area index (also named as the leaf area density) in the vertical direction of the ccontinuous canopy, the equation for which is: Here, L is the leaf area index (LAI) and f (θ l ) is the leaf inclination distribution function (LADF), which uses an ellipse distribution [51] (This is an equation for azimuthally independent analysis [19]).
Since E s , E − , E + , and E o are assumed to be vertical in the four-stream radiative transfer equations (Z-axis in Figure 1d), Equation (1) describes the extinction of specular flux in the vertical direction, Equations (2) and (3) describe the vertical transfer processes of the upward and downward diffuse fluxes, and Equation (4) is an approximation of the radiative transfer equation, referring to the entire vertical radiative transfer process from the viewing [52]. Combined with the assumption that the canopy is infinitely extended (Figure 1c) for four-stream radiative transfer equations, radiative transfer of vertical flux density achieves good results in continuous crops. However, for the row canopy, there is a radiative transfer parallel to the horizontal ground surface and passing through the lateral "walls" (Figure 1a). Therefore, the lack of horizontal radiative transfer may cause an energy imbalance when crops have a row structure, which affects the accuracy of the model.

Refined Four-Stream Radiative Transfer Equations for Row-Planted Crops Considering Horizontal Radiative Transfer Equations
The radiation exchanged between the lateral "walls" and the bottom of the atmospheric boundary layer is defined as horizontal radiation, which is the diffuse radiation parallel to the horizontal ground surface and passing through the lateral "walls". There were two lateral "walls" in the canopy closures. For one type of lateral "wall", the emitted fluxes included specular flux and diffuse flux, which we Remote Sens. 2020, 12, 1290 5 of 28 named lateral "wall" A ("bc" Figure 1a). For the other type of lateral "wall", the emitted flux was only the diffuse flux, which we named lateral "wall" B ("ad" and "ef" in Figure 1a). To improve the radiation energy balance issue in row-planted crops, we needed to consider this horizontal radiative transfer of the lateral "walls" into the original four-stream radiation transfer equations. According to the derivation in Supplementary Material B, we derived the horizontal radiative transfer equations, i.e., Equations (B-12), and (B-13). Combining horizontal radiative transfer equations with Equations (1)-(4), we refined the four-stream radiative transfer equation into: The optical coefficients (i.e., k, K, s , s, σ, w, v, v , and a) in Equations (6)-(9) used the previous equations in [19]. For the optical coefficients in the horizontal radiative transfer equations (Equations (10) and (11), the equations are: Here, n is the attenuation coefficient for horizontal diffuse flux, g is the radiative converted coefficient describing the proportion of downward diffuse flux converted to horizontal diffuse flux from the lateral "wall"; g is the radiative converted coefficient describing the proportion of upward diffuse flux converted to horizontal diffuse flux of the lateral "wall"; m is the bidirectional scattering coefficient from specular flux to horizontal diffuse flux (for their derivation information see Supplementary Material B). For unknowns in Equations (10) and (11), E b is diffuse horizontal hemispheric flux density through lateral "wall" A and E d is diffuse horizontal hemispheric flux density through lateral "wall" B. Finally, L row is the differential leaf area index for canopy closure in the vertical direction, U is the horizontal differential leaf area index, and their equations are: Here, A 1 is the row width, A 2 is the distance between rows, and h is the height of the canopy. According to the definition of E b and E d , E b and E d are horizontal radiation in lateral "wall" A and Remote Sens. 2020, 12, 1290 6 of 28 lateral "wall" B, respectively. Therefore, Equation (10) is the horizontal radiative transfer equation in lateral "wall" A and Equation (11) is the horizontal radiative transfer equation in lateral "wall" B. These two equations are considered, and the entire system of equations satisfies the radiation energy balance for the canopy closure of row-planted crops.
We used z = 0 and z = −1 as the boundary conditions in the vertical direction, and B = A 1 2h sin ϕ r − ϕ o (it's acquired in Supplementary Material B-2) as the boundary conditions in the horizontal direction to solve Equations (6)- (11). We calculated the directional reflectance factor (DRF) on the boundary of the canopy closure; that is, the DRF at the top of the canopy closure (R c ) and the DRFs of lateral "wall" A (R b ) and lateral "wall" B (R d ). Among them, for the DRF at the top of canopy closure (R c ), we referred to the solution in the original four-stream radiative transfer equations (i.e., Equations (29)-(33) in Section 2.3.1). For the DRFs of lateral "wall" A (R d ) and lateral "wall" B (R d ), we showed the derivation in Supplementary Material D. The equation for the DRF of lateral "wall" A is: (18) and the equation for the DRF of lateral "wall" B is: Here, E s (B) is the specular flux on the surface of the lateral "wall" and E ± (B) is the diffuse flux on the surface of the lateral "wall". Their equations are: Here, r sd , r do , r so , and r dd in Equations (18) and (19) are: With the calculated results of these DRFs, we then combined the area fractions of each component in row crops calculated from the gap probabilities (Supplementary material C) and DRF at the top of between row area (R br ) (it was calculated from the integral radiative transfer equation in Supplementary material E) to build a row-planted crop model to calculate the DRF.

Sum of the Directional Reflectance Factor of Row-Planted Crops
We followed the assumption that the canopy extended infinitely in the four-stream radiative transfer equations [19], and the row canopy also alternated between canopy closure and between row to reach infinite extension. Therefore, the DRF at the top of canopy closure (R c ) and the DRF at the top of the between row area (R br ) composed the sum DRF of row-planted crops, namely: According to the composition of the DRFs in Equation (26), we explained the calculation for DRF at the top of the canopy closure and at the top of the between row area, respectively.

Directional Reflectance Factor of the Canopy Closure
In the field of radiative transfer, the DRF at the top of the canopy closure included the single scattering (R c_1 ) and the multiple scattering of the canopy closure (R c_m ), the equation for which is: Below is the DRF derived from the original four-stream radiative transfer equations in [13]: This equation shows that the DRF includes the bidirectional reflectance factor (r * so ) and the hemispherical directional reflectance factor (r * do ), which is the ratio of the reflected flux (r * so E s (0) + r * do E − (0)) to the incident flux (E s (0) + E − (0)) at the top of the canopy. The canopy closures are different from the continuous canopy-they are separated by between rows (Figure 1a). In this section, the differential leaf area index (leaf area density) (L row ) and the area fractions of each component (S) changed according to the row structure. Therefore, to solve the canopy closure using the solution method in the original four-stream radiative transfer theory, L row and S needed to be modified according to the structure. The differential leaf area index (leaf area density) for canopy closure in the vertical direction was modified in Equation (16), and the area fractions of each component were modified in Table C Single Scattering of the Canopy Closure According to [13], we used the following equation to calculate the single scattering of the canopy closure (R c_1 ): (29) in which: Here, r 1 so_v is the single scattering of specular flux in the canopy closure, r 1 so_s is the single scattering of specular flux from the soil in the canopy closure, and r s is the DRF of the soil. S closure_s (z) is the fraction of observed canopy illuminated by the specular flux in the canopy closure, S closure_s (h) is the fraction of observed soil illuminated by the specular flux in the canopy closure, and their expressions are shown in Table C Multiple Scattering of the Canopy Closure According to [13], we used the following equation to calculate the multiple scattering of the canopy closure: (31) in which: Here, r m so_v is the multiple scattering of specular flux in the canopy closure, r m so_s is the multiple scattering of specular flux between soil and vegetation in the canopy closure, r 1 do is the single scattering of diffuse flux in the canopy closure, and r m do is the multiple scattering of diffuse flux in the canopy closure. S closure_d is the fraction of canopy closure illuminated by the diffuse flux, and its expressions are shown in Table C-1 in Supplementary Material C.

Directional Reflectance Factor of the Between Row Area
Similar to Section 2.3.1, the DRF at the top of the between row area included the single scattering of the between row area (R br_1 ) and the multiple scattering of the between row area (R br_m ), the equation for which is: Single Scattering of the Between Row Area The single scattering of the between row area is: (35) in which: Here, S between_row_s is the fraction of observed soil background in the between row area illuminated by the specular flux, and its expressions are shown in Table C-1 in Supplementary Material C.

Multiple Scattering of the Between Row Area
The multiple scattering of the between row area is: Here, R s_m is the multiple-scattering of soil in the between row area, which is: We calculated this parameter based on the solution for the integral form radiative transfer equation. For detailed calculation, please refer to Supplementary Material E. S between_row_d is the fraction of the between row background illuminated by the diffuse flux. Its expressions are shown in Table C-1 Remote Sens. 2020, 12, 1290 9 of 28 in Supplementary Material C. We used the horizontal DRFs (i.e., R b and R d ) as the initial values to calculate the multiple scattering of the between row area.

Input Parameters of the MFS Model
Based on the above mathematical analysis, we refined the four-stream radiative transfer equations by considering the horizontal radiative transfer equations, and then solved them to build the model for row-planted crops combined with the reflectance issue of the between row reflectance, then proposed the MFS model. In the MFS model, the directional hemisphere reflectance and transmittance of leaves are provided by the PROSPECT-5 simulations [53,54]. The input parameters of MFS include: (1) Geometrical parameters: solar zenith angle (θ s ), solar azimuth angle (ϕ s ), viewing zenith angle (θ o ), viewing azimuth angle (ϕ o ), and row azimuth angle (ϕ r ); (2) Canopy parameters: height of the canopy (h), row width (A 1 ), distance between rows (A 2 ), leaf area index (L) and effective leaf area index (L E ), average leaf inclined angle (θ l ), and canopy dimension parameter (l * L ); (3) Biochemical leaf parameters: chlorophyll content (C ab ), carotenoid content (C ar ), brown pigment content (C brown ), equivalent water thickness (C w ), leaf mass per unit leaf area (C m ), and structure coefficient (N); (4) Canopy radiative parameters: the fraction of incoming diffuse radiation (skyl).

Data
In this study, we validated the MFS model using both computer simulation and in situ measurements.

Computer-Simulated Data
In the study of the canopy radiative transfer model, the simulation of the canopy radiative transfer assumed that there was no influence from the atmosphere. Since the in situ measurement is done in the atmospheric environment, it is difficult to exclude the influence of essential climate variables (ECVs) on the "true value" of the measurement [55]. According to the radiation transfer model intercomparison (RAMI), the DRF calculated by 3D computer simulation has the highest accuracy, and the calculation result can be used as a "surrogate truth" to validate the turbid medium model in radiative transfer [56][57][58]. The MFS model in this study is based on the radiative transfer equations of flux density in order to calculate the reflectance. Therefore, we used a 3D computer simulation based on the principle of flux density to validate the MFS model. The 3D computer simulation is the extended a 3D radiosity graphics (RGM) model [11,59]. The RGM model is a radiosity model based on a bilinear equation (a simple non-linear differential equation) [5]. It uses numerical calculation methods to calculate the scattering of polygons in a scene constructed by a computer graphics method and has high calculation accuracy [11]. To use the RGM model to validate the MFS model, we first needed to use computer graphics to construct the abstract scenes that were assumed to have medium turbidity, and then set the inputs and outputs for the two models.

Generation of Computer Abstract Scenes
According to the study in the fourth phase of the radiative transfermodel intercomparison (RAMI-IV, http://rami-benchmark.jrc.ec.europa.eu), the establishment of the scene is based on measurements [60]. Therefore, measurements of corn during the growth period in the Yingke Oasis (see Section 3.2.1 for details) were used to establish the abstract scenes [61]. The abstract scene was generated (Figure 2), in which the canopies with four geometric structures were generated to reflect the row growth state, including the proportion of between row dominance (stage_rv1), proportions of between row and canopy closure equality (stage_rv2), the proportion of canopy closure dominance (stage_rv3), and continuous vegetation (stage_cv). The parameters of the constructed abstract scene are listed in Table 1. validate the MFS model, we first needed to use computer graphics to construct the abstract scenes that were assumed to have medium turbidity, and then set the inputs and outputs for the two models.

Generation of Computer Abstract Scenes
According to the study in RAMI-IV (http://rami-benchmark.jrc.ec.europa.eu), the establishment of the scene is based on measurements [60]. Therefore, measurements of corn during the growth period in the Yingke Oasis (see Section 3.2.1 for details) were used to establish the abstract scenes [61]. The abstract scene was generated ( Figure 2), in which the canopies with four geometric structures were generated to reflect the row growth state, including the proportion of between row dominance (stage_rv1), proportions of between row and canopy closure equality (stage_rv2), the proportion of canopy closure dominance (stage_rv3), and continuous vegetation (stage_cv). The parameters of the constructed abstract scene are listed in Table 1.  Table 1. Values of parameters in the constructed abstract scene. Here, the shape of the leaf is triangular, with the area ranging from 0.0002 m 2 to 00003m 2 in the scenes; the shape of the soil area is square, with an area of 0.0025 m 2 ; and the leaf inclined angle is randomly distributed in the zenith and azimuth directions. The symbols in Table 1 are: A 1 = row width, A 2 = distance of between-rows, h = the height of the canopy, L = leaf area index, n = the number of triangular leaves, w* = one of the short sides of the horizontal triangle leaf, l* = the other short sides of the horizontal triangle leaf, θ l is average leaf inclined angle. ϕ r is row azimuth angle.

Input and Output Settings for Computer-Simulated Validation
Input Setting In both the RGM model and MFS model, the directional hemisphere reflectance and transmittance of leaves were provided by the PROSPECT-5 model. The input parameters of PROSPECT-5 model (C ab = 40.29 ug·cm −2 , C ar = 8.54 ug·cm −2 , C brown = 0, N = 1.52, Cw = 0.014 cm −2 , C m = 0.004 g·cm −2 ) were also obtained (see Section 3.2.1 for details). The soil reflectance was obtained by measurement, in which the red band with 670 nm was 0.242 and the near-infrared band with 850 nm was 0.314. According to data from the Yingke Oasis, the sun zenith angle was 25 • and the sun azimuth angle was 130 • . To keep the flux density in both the MFS model and RGM model consistent, we changed the measured flux density to the normalized flux density in the validation. In the MFS model, this took into account the diffuse incidence (Equation (33)) at the top of the canopy that was ignored by the previous row model [40], for which the fraction of incoming diffuse radiation (skyl) for both models was set to 0.1. Combined with Equation (C-22) in Supplementary Material C, n ∆ , l * , and w * in Table 1 were used to calculate the canopy dimension parameter (l * L ) in the MFS model. Finally, for the row structure and in the MFS model, we used the measured values of the Yingke Oasis (i.e., the values of the abstract scene in Table 1).

Output Setting
Previous studies showed that there was a systematic error in the four-stream radiative transfer equations when the viewing azimuth angle was greater than 40 • , and this phenomenon was particularly obvious in the NIR band [40]. To better show that the MFS model had solved this problem, we showed the DRF field of the entire hemisphere space-simulations were performed with the viewing azimuth angle varying from 0 to 350 • with a step of 10 • , and with the viewing zenith angle ranging from 0 to 89 • with a step of 5 • . One thing needs to be explained-the MFS model assumed that the canopy was an infinite extension. To unify the comparison conditions, we chose the infinite canopy mode in RGM to calculate the DRF. The infinite canopy does not just calculate the DRF of the finite canopy for the 1-to 2-row cycle shown in Figure 2, but copies the finite canopy for 101 × 101 times as a duplicated canopy to obtain an approximately infinite extension, and then calculates the DRF of this infinite canopy (shown in Appendix B in [11]) as its calculation principle. When the viewing zenith angle was 90 • , the RGM model was unable to be simulated. Therefore, the maximum angle is 89 • . Finally, Figures 3,5,7,and 9 show the DRF fields simulated by both models, while Figures 4,6,8,and 10 show the corresponding four viewing modes for row-planted crops in DRF fields, i.e., the principal plane (PP, viewing azimuth angle is aligned with solar azimuth angle), orthogonal plane (OP, viewing azimuth angle is perpendicular to solar azimuth angle), along row plane (AR, the plane is along the row direction), and orthogonal row plane (OR, the plane is perpendicular to the row direction).

Measurement Experiment
The in situ data were collected in the Watershed Allied Telemetry Experiment Research (WATER) project, ranging from May 20 to July 9, 2008, in the site located at 38.8571 • N, 100.4104 • E, at an elevation of 1529 m [62,63]. This is a cornfield located in the Yingke Oasis, Zhangye City, Gansu Province, China. The area of the plots is 180 × 180 m. In each plot, three samples were randomly selected to perform measurements [64]. The measurements included two of the crop canopy, one of row-planted crops at the jointing stage, and the other of continuous crops at the heading stage. The data used in this study were introduced as outlined below.

Directional Reflectance Factor (DRF)
The DRF of the corn and the DRF of soil at wavelengths ranging between 400 and 2500 nm were measured using the Spectrometer (It name is ASD FieldSpec Pro Spectrometer) mounted on a multi-angle observation frame. The distance from the sensor to the bottom of the canopy was 5 m, and the viewing angle under FOV was 25 • during the measurement. The viewing azimuth angle varied from 0 • to 360 • along the azimuth of the plane (i.e., four viewing modes described below), while the zenith angle varied from −60 • (i.e., backward direction) to 60 • (i.e., forward direction) with an interval of 10 • . For further analysis of the anisotropy of DRFs for row-planted crops, these measurements were separated into four categories according to four viewing modes, namely PP, OP, AR, and OR.

Canopy Structural and Leaf Biochemical Parameters
The effective leaf area index and average leaf angle were measured using the LAI-2000. The leaf area index, row width, distance between rows, average canopy height, row azimuth angle, average width of leaves, average length of leaves, and average area of leaves were simultaneously measured using the direct measuring method (the leaf area index was measured using the improved harvesting method-IHM [65]) [61]. The chlorophyll content was calculated using the empirical formula based on the SPAD index (this is a vegetation index that determines the chlorophyll in the leaf), acquired using the SPAD chlorophyll meter [66,67]. The leaf mass per unit leaf area (C m ) and the equivalent water thickness (C w ) were derived from the weights of fresh and dry leaves [68], given in g·cm −2 [69]. The carotenoid content (C ar ), brown pigment content (C brown ), and structure coefficient (N) were determined using the Lopex1993 database [70].  Table 2 showed other input parameters in the three models. Here, the abbreviations in the title are: MFS = modified four-stream radiative transfer model; DRM = a spectral directional reflectance model; SAIL = light scattering by leaf layers with application to canopy reflectance modeling. The symbols in Table 2

Output Setting
To validate whether the MFS model improves the accuracy of canopy reflectance of row-planted crops at zenith angles greater than 40 • , we chose θ o = ±60 • and θ o = ±50 • as quasi-horizontal DRFs to be used in the validation. On this basis, to further reflect the reflectance anisotropy issue, we selected three measurement modes that lateral "walls" may observe, namely PP, OP, and OR. To further show that the MFS model had improved the accuracy of the canopy reflectance when the viewing azimuth angle was greater than 40 • compared with the previous model, we chose the original row model for comparison, i.e., the DRM model [40]. To further assess the applicability of the model for the continuous crops, simulations of the SAIL model were also used for comparison with the simulations of the MFS model in the vertical viewing direction (θ o = 0 • ). Finally, Figures 11 and 12 show the corresponding DRF simulations and in situ measurements.  Figure 3). The contours of the fields are slightly different. However, the magnitude of the DRF fields simulated by the two models is between 10 −3 and 10 −2 , and the difference between the two model simulations is not more than 16.69% (Table 3). Hence, the difference is very small. From Figure 5a,b to Figure 5g,h, the high-value area generated by the soil background shrinks into hotspots as the crop grows. In the cross-sectional view of the DRF fields, i.e., Figure 4, the difference in the contours of the fields in Figure 3 is very small. For the four viewing modes in the red band, DRFs simulated by the MFS model and DRFs simulated by the RGM model are more consistent, except for the simulation at a viewing zenith angle greater than 80 • along row plane (AR) mode (Figure 4c).

Validation of MFS Model Using Computer-Simulated Data
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 33     Figure 5), and the difference between the two models does not exceed 11.69% (Table 3). Different from the shape of DRF fields in the red band, the shape of DRF fields in the NIR band presents a bowl-shape ( Figure 5). In the growth stage of the DRF fields, the bowl edges (the area where the reflectance increases rapidly) shrink toward a small zenith angle, e.g., the bowl edge is around 80° for stage_rv1 ( Figure 5(a-b)), around 70° for stage_rv2 ( Figure 5(c-d)), around 60° for stage_rv3 ( Figure 5(e-f)). From the AR direction, i.e., stripe area from 0 to 180 ° in the DRF fields, there is a distinct low-value area. Among the cross-sectional views of the four viewing modes in Figure 5, namely Figure 6 Figure 5), and the difference between the two models does not exceed 11.69% (Table 3). Different from the shape of DRF fields in the red band, the shape of DRF fields in the NIR band presents a bowl-shape ( Figure 5). In the growth stage of the DRF fields, the bowl edges (the area where the reflectance increases rapidly) shrink toward a small zenith angle, e.g., the bowl edge is around 80 • for stage_rv1 (Figure 5a,b), around 70 • for stage_rv2 (Figure 5c,d), around 60 • for stage_rv3 (Figure 5e,f). From the AR direction, i.e., stripe area from 0 to 180 • in the DRF fields, there is a distinct low-value area. Among the cross-sectional views of the four viewing modes in Figure 5, namely Figure 6  Unlike the red band, the multiple scattering effects of the NIR band can be ignored. For the NIR band, it is very important. To further analyze the differences between the RGM model and the MFS model in the simulations, DRFs in the NIR band are divided into two groups, i.e., single scattering and multiple scattering DRFs. For single scattering, the difference between the two models does not exceed 9.99% (Table 3). When the zenith angle is greater than 40° (outside of shaded annulus in Figure 7), the consistency between the DRFs simulated by the two models is very high (Figure 7). When the row structure is obvious (Figure 7(c-f)), the low-value area for the DRFs appears in the opposite direction to the incident hemisphere (   Figure 8a were slightly different, the overall consistency was high. Figure 9 displays changes in the growing stage of the crops due to the multiple scattering of the NIR band. This shows that the DRFs simulated by the MFS model are slightly higher than DRFs simulated by the RGM model, and the difference between the two models does not exceed 3.35%, except for stage_rv1 (Figure 9a,b). With changes in the row structure, the low-value area shrinks toward the center of the row angle. For the cross-sections of the four viewing modes in Figure 10  Unlike the red band, the multiple scattering effects of the NIR band can be ignored. For the NIR band, it is very important. To further analyze the differences between the RGM model and the MFS model in the simulations, DRFs in the NIR band are divided into two groups, i.e., single scattering and multiple scattering DRFs. For single scattering, the difference between the two models does not exceed 9.99% (Table 3). When the zenith angle is greater than 40 • (outside of shaded annulus in Figure 7), the consistency between the DRFs simulated by the two models is very high (Figure 7). When the row structure is obvious (Figure 7c-     Here VZA = viewing zenith angle. Figure 9 displays changes in the growing stage of the crops due to the multiple scattering of the NIR band. This shows that the DRFs simulated by the MFS model are slightly higher than DRFs simulated by the RGM model, and the difference between the two models does not exceed 3.35%, except for stage_rv1 (Figure 9a,b). With changes in the row structure, the low-value area shrinks toward the center of the row angle. For the cross-sections of the four viewing modes in Figure 10, the DRFs of the multiple scattering simulated by the MFS model are highly consistent with the DRFs of the multiple scattering simulated by the RGM model.  The white shaded annulus in Figure 9. (g,h) is used to distinguish it from the DRF field.
Remote Sens. 2020, 12, x FOR PEER REVIEW 24 of 33 the north direction, the blue line represents the PP mode, and the shaded annulus lies at an azimuth angle of 40°. (g,h) The white shaded annulus is used to distinguish it from the DRF field.

Quantitative Analysis Results
The results of the qualitative analysis between the two models are shown in Section 4.1.1 and Section 4.1.2. To further describe the comparison between the two models, we present the results of the quantitative analysis. Figures 3, 5, 7, and 9 include 684 samples simulated under the same sun and viewing geometries for each abstract scene. The average difference of the DRF fields between the two models is summarized in Table 3. In Table 3, the average difference between the DRF simulated by the MFS model and the DRF simulated by the RGM model is generally not more than 10%, except for DRF_red (-16.69%) and DRF_NIR_m (-22.64%) in the stage_rv1. On the same computer with the same configuration, the calculation time for the MFS model was 0.135 hours for all of the scenes, while the calculation time for the RGM model was 5.678 hours for all of the scenes. In Table 3, the DRFs simulated by the MFS model show good correlation with DRFs simulated by the RGM model, with correlation coefficients greater than 0.9341, except for the DRFs in the NIR band in stage_cv (R=0.8207). RMSEs are less than 0.0019 in the red band and less than 0.0187 in the NIR band. The MFS model simulates the DRFs contributed by single scattering and multiple scattering in the NIR band with good accuracy, i.e., R is greater than 0.96 and RMSEs are less than 0.0088 for single scattering, while R is greater than 0.9 and RMSEs are less than 0.012 for multiple scattering.

Quantitative Analysis Results
The results of the qualitative analysis between the two models are shown in Sections 4.1.1 and 4.1.2. To further describe the comparison between the two models, we present the results of the quantitative analysis. Figures 3, 5, 7 and 9 include 684 samples simulated under the same sun and viewing geometries for each abstract scene. The average difference of the DRF fields between the two models is summarized in Table 3. In Table 3, the average difference between the DRF simulated by the MFS model and the DRF simulated by the RGM model is generally not more than 10%, except for DRF in red band (DRF_red) (−16.69%) and multiple scattering of DRF in NIR band (DRF_NIR_m) (−22.64%) in the stage_rv1. On the same computer with the same configuration, the calculation time for the MFS model was 0.135 hours for all of the scenes, while the calculation time for the RGM model was 5.678 hours for all of the scenes. In Table 3, the DRFs simulated by the MFS model show good correlation with DRFs simulated by the RGM model, with correlation coefficients greater than 0.9341, except for the DRFs in the NIR band in stage_cv (R = 0.8207). RMSEs are less than 0.0019 in the red band and less than 0.0187 in the NIR band. The MFS model simulates the DRFs contributed by single scattering and multiple scattering in the NIR band with good accuracy, i.e., R is greater than 0.96 and RMSEs are less than 0.0088 for single scattering, while R is greater than 0.9 and RMSEs are less than 0.012 for multiple scattering.

Validation of the MFS Model Using In Situ Data
Figures 11a-f and 12a-f show the simulation differences between the row models under the PP, OR, and OP at zenith angles of 60 • and 50 • , respectively. MFS 1 and MFS 2 are highly consistent with the measured data. The DRFs simulated by MFS 2 are better than the DRFs simulated by MFS 1 , and their correlation coefficients (R) are higher than 0.985 with RMSEs of less than 0.023 ( Table 4). The DRF simulated by the DRM model is seriously overestimated when the zenith angle is greater than 40 • , especially in the NIR band (Figure 11a-f). Comparatively, the DRF simulated by MFS 1 and MFS 2 at the zenith angle of 60 • is slightly better than that at the zenith angle of 50 • . For the continuous crops, the performance of MFS and SAIL models is comparable, and DRFs simulated by the two models have the same R and RMSE values (Figure 11g).     Here RC = row-planted crops, CC = continuous crops, PP = principal plane, OP = orthogonal plane, OR = orthogonal row.

Systematic Deviation of the DRFs in the NIR Band at Zenith Angles Larger than 40 •
Comparing the validation using computer simulation and the validation using in situ measurements, the MFS model had higher accuracy in calculating DRFs. For the validation using computer simulation, the MFS model achieved high accuracy during the growth stage of crops (Figures 3, 5, 7 and 9) and effectively reduced the systematic deviation of the DRFs in the NIR band calculated by the row model based on the four-stream radiative transfer theory when the zenith angle was greater than 40 • (outside of shaded annulus in Figure 5). For validation using in situ measurements, the consistency of the MFS model (MFS 1 and MFS 2 ) and the in situ measurement was maintained to the greatest extent and solved the systematic deviation of the DRM model in the NIR band when the zenith angle was greater than 40 • (Figure 11a-f). Compared with MFS 1 and MFS 2 , the clumping index and the canopy dimension parameter (l * L ) did not obviously improve the overall accuracy (Figures 11a-f and 12a-f). Therefore, when the zenith is larger than 40 • , the large systematic deviation of the DRF in the NIR band simulated by the DRM model is not caused by these two factors. A previous study on the DRM model pointed out that the systematic deviation of the DRFs in the NIR band at zenith angles larger than 40 • may be caused by the coefficient calculation method SAIL model used by the four-stream radiative transfer equations [40]. In this study, the calculation method for the coefficients (k, K, s , s, σ, w, v, v , and a in Equations (1)-(4) and Equations (6)-(9), i.e., SAIL model [19]) was used to solve the modified four-stream radiative transfer equations, but the systematic deviation of the DRFs in the NIR band did not occur when the zenith angle was larger than 40 • (Figures 11a-f and 12a-f). These results also showed that when the zenith angle is greater than 40 • , the large systematic deviation of the DRFs in the NIR simulated by the DRM model comes from other aspects. We further analyzed the core equations in the DRM model, but the horizontal radiative transfer was not considered. For the equation of the canopy closure for multiple scattering in the DRM model (i.e., Equation (13) in [40]), this was derived from an approximate radiative transfer equation [52] (i.e., Equation (4) in this study; the derivation process is detailed in Appendix A in [40]) instead of all of the four-stream radiative transfer equations (i.e., Equations (1)-(4)). These mathematical derivations imply that the original four-stream radiative transfer equations are subjected to a second mathematical approximation. Therefore, the calculation accuracy of the DRFs is reduced in the multiple scattering of the canopy closure. Based on the above analysis, the systematic deviation of the DRFs in the DRM model (Figure 11a-f) may come from the multiple scattering equation for the canopy closure and the lack of horizontal radiative transfer. The MFS model did not make any assumptions or simplifications in the original four-stream radiative transfer equations, directly calculated the exact solutions for the differential equations, solved the MFS radiative transfer equations (Section 2.2 and Supplementary material D), and then considered the DRF for between rows (Supplementary materials E) to build the MFS model, in which it was applied in the multi-angle simulation of DRF fields in row-planted crops.

Differences between MFS Model and Computer Simulation
The DRF values simulated by the MFS model produced results that were essentially in agreement with the DRF values simulated by the RGM model. However, the contours in the DRF fields simulated by the MFS model did not have the jitter of the contours in the DRF field simulated by the RGM model ( Figures 3, 5, 7 and 9). The same phenomenon also occurred in [40]. For the differences in the contours in the DRF fields, these should come from the systematic deviation of the algorithm between the two models. Further analyzing the algorithms, in the MFS model, the amplitude of the DRF fields should be determined by the area fraction of each component calculated by the gap probabilities (Supplementary Material C) and the solution method for the modified four-stream radiative transfer equations (non-numerical adding method, i.e., R c , R b , and R d in Supplementary Material D and R br in Supplementary Material E). In the RGM model, the amplitude of the DRF fields should be determined by the visible polygons factors (i.e., visual probabilities of two polygon elements) [11] and the numerical algorithm (Gauss-Saidel method) to solve the system of equations [59]. As a result, the differences in results are generated based on the differences in the core solutions described above.
In further analysis of the four anisotropic DRF viewing modes (PP, OP, AR, and OP), DRFs simulated by the MFS model and DRFs simulated by the RGM model were more consistent, except when the viewing zenith angle was greater than 80 • in the AR mode (Figure 4c). The reason for this phenomenon may be the difference between the design strategies for the extinction process in the two models. To calculate extinction, the MFS model used a penetration function based on optical paths (Equation (C-5) in Supplementary Material C), and RGM used the visual factors between polygons and used a duplicated canopy (Section 3.1.2). The calculation difference between the two extinction principles was the main cause of the system deviation. When the fraction of the distances between rows was large, the systematic deviation occurred between the two models, e.g., the maximum absolute value for the average difference in stage_rv1 was 22.64% (Table 1), while the maximum absolute value for the average difference in stage_rv2 was 13.76% (Table 1). The possible reason is that the two models use different algorithms to calculate the DRF of the soil between rows. The MFS uses a one-dimensional integral equation (Equation (E-1) in Supplementary Material E), which calculates the isotropic DRF of the soil between rows. The RGM is a three-dimensional computer model that calculates the anisotropic DRF of the soil between rows (e.g., the brown squares shown in Figure 2 are approximated as an anisotropic soil). Therefore, as the distances between rows increases, the area fraction of soil increases, and the differences in DRFs become more obvious.
Based on the analysis of these algorithms, aside the RGM model needing to use computer graphics technology to construct computer scenes, the biggest difference between the MFS model and the RGM model was the algorithms used to solve the equations. The MFS model used a non-numerical method to solve the radiative transfer equation, while the RGM model used a numerical method to solve the radiative transfer equation. The accuracy of the numerical algorithm was very high, but achieving this higher accuracy took more computing time (e.g., 5.678 hours for the RGM model, shown in Section 4.1.3). Therefore, computer simulation is mostly used for the understanding of the radiation regime and as a "surrogate truth" to validate other models [5]. Although the radiative transfer model used a non-numerical method with a small loss in calculation accuracy, the calculation time was greatly reduced (e.g., 0.135 hours for the MFS model, shown in Section 4.1.3), making it the best choice for the ill-conditioned inversion issue, which required large datasets.

The Relationship between the MFS Model and SAIL Model
As an extension of the four-stream radiative transfer equations in row-planted crops, the horizontal radiative transfer (Equations (10) and (11)) was considered in the MFS model, and the non-numerical method for four differential equations with four-stream radiative transfer theory continued to be used. For continuous crops, the performance of the MFS and SAIL models is comparable, while DRFs simulated by the two models had the same R and RMSE values (Figure 11g and Table 4). The reason for this is that the MFS model used a form of overdetermined systems (six equations with only four unknowns) in mathematics, the horizontal radiative transfer equations were two redundant equations, and these were degraded in the simulation for continuous crops, which seamlessly connected the MFS radiative transfer equations and original four-stream radiative transfer equations using mathematics. This implied that the MFS model extended the original four-stream radiative transfer equations in row-planted crops.

Conclusions
In this study, we presented a refined radiative transfer model for row-planted crops. We proposed a modified four-stream (MFS) radiative transfer model by considering horizontal radiation through the lateral "walls" and took row structure into consideration for more accurate estimation. We validated the MFS model using both computer simulations and in situ measurements. This shows that the MFS model can be used to simulate crop canopy reflectance at different growth stages, with an accuracy comparable to the computer simulations. Moreover, the MFS model can be successfully used to simulate the reflectance of continuous and row crop canopies, and therefore can address the large viewing zenith angle problems in the previous row model based on four-stream radiative transfer equations, especially when the zenith angle is greater than 40 • . Our results demonstrate that horizontal radiation is an important factor that needs to be considered in modeling the canopy reflectance of row-planted crops. Hence, the refined four-stream radiative transfer model is applicable to the real world.
Author Contributions: X.M., T.W., and L.L. conceived the idea. X.M. completed the mathematical derivation, data analysis, and drafted the manuscript. All authors contributed to the editing of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
The in situ data used in this paper were supported by the Watershed Allied Telemetry Experiment Research (WATER) project. This work was supported by the National Natural Science Foundation of China (Grant 41401372).