**4. Discussion**

### *4.1. Multiple Scattering of Row Crops in the GO Approach*

This study substantiates similar findings obtained in earlier studies by showing that it is necessary to consider the multiple-scattering contribution in the NIR band in the GO approach [13,26,27]. Moreover, we considered the multiple-scattering contribution in row crops as an extension of the

previous study field (forest field) to the object of row crops. In the row modeling based on the GO approach, its initial purpose was to provide a physical basis to calculate temperature (or emissivity). However, due to the different physical properties between temperature (or emissivity) and reflectance, the previous studies rarely considered multiple scattering issues in row modeling [14,28,29]. Although subsequent studies attempted to consider multiple scattering issues in the reflectance modeling of row crops, these models usually consider a row structure to modify existing continuous crop models based on the RT approach [17,31,54,55]. Therefore, their physical basis was still the volume scattering of the radiative transfer equation, and their essence was the RT model. Different from the previous modeling, our model took the multiple-scattering equation constructed by the RT approach in GO modeling further to modify the GO approach. The multiple-scattering equation achieved a good calculation accuracy (Figure 7b,d,f,h). In the multiple-scattering contributions of row crops, the radiation energy of multiple scattering in the NIR band was very strong, accounting for a large proportion of the sum of the reflectance (comparing Figure 6b,d,f,h and Figure 7b,d,f,h). Specifically, as crops gradually grew from row crops to continuous crops, the proportion of multiple scattering in the NIR band in the sum of the reflectance gradually increased from 20% to 50% (comparing Figure 6b,d,f,h and Figure 7b,d,f,h). These results imply that the sum of the reflectance would be significantly underestimated if the multiple-scattering contribution is not considered, and the underestimation becomes more obvious as LAI gets bigger. Therefore, the multiple-scattering equation is considered in GO modeling and can be accurately calculated, which can effectively improve the accuracy of the sum of the reflectance in the NIR band. However, in Table 4, the mean absolute deviation (MAD) of R\_NIR\_m gradually decreased from 9.82% in the Stage\_rc1 to 0.99% in the Stage\_cc. This result implies that the row structure influences the calculation accuracy of multiple scattering, whereby, when the row width decreases, the between-row distance increases, and the difference in the multiple-scattering contribution between our model and the RGM model will increase. This conclusion is consistent with the simulation results in [30]. The primary reason is that our model describes Lambertian soil (soil reflectance as a single value, see Equations (13) and (22)), but the RGM model describes non-Lambertian soil (the brown squares shown in Figure 4 are approximated as an anisotropic soil, hence soil reflectance is a multivalue consisting of multiple squares). Therefore, the increase in the between-row distance will cause a difference in the multiple scattering between our model and the RGM model.

In the GO approach, we studied an adding method (Section 2.1.2) and a mathematical solution of integral radiative transfer equation (Section 2.1.3) to derive the multiple-scattering equation. This is a typical method for the multiple-scattering equation, which was established by the RT approach, to be coupled into the GO approach for a modified GO model. This method is based on the RT approach to derive the multiple-scattering equation. However, in a previous study (mainly in the heterogeneous forest), this method often ignored the geometric characteristics between individual canopies [26], which was reflected in the results as a computational deviation in the sum of the reflectance [13,27]. The between-individual canopies mentioned in [13,27] are the between-row area in row crops. Our research focused on this point in row modeling. In our model, the openness angle of the between-row area (α1) and nonopenness angle of the between-row area (α2) were introduced to calculate an average value of the escape probability of radiation in the between-row area (*P*open). The above treatment showed the scattering direction (or scattering angle with typical geometric characteristics) in the multiple-scattering contribution of the between-row area, and further reflected the influence of geometric characteristics (*A*2, ϕ*r* and *h* in Equation (20)) on the multiple-scattering equation (Equation (22)). The multiple-scattering contribution simulated by our model had higher R-values and lower RMSEs compared to the RGM model (Table 4). The modified treatments of these equations and results further illustrated that our model considered how the geometric characteristics between individual canopies can improve the calculation accuracy of the multiple-scattering contribution. We further analyzed the sum of the reflectance simulated by our model compared to those of the RGM model was of high consistency (Figure 6). The simulation results and in situ measurement of our model also have good accuracy overall (Figures 9–11). These results show that geometric characteristics

considered in the multiple-scattering equation established by the RT approach are necessary for GO modeling. This consideration not only improves the accuracy of multiple scattering but also the accuracy of the sum of the reflectance.

When we analyzed the multiple-scattering equation in the previous model in the heterogeneity forest (Equation (27) in [27] and Equations (8)–(11) in [13]), we found that the single-scattering albedo of the leaf was added to derive the multiple-scattering equation. The purpose of this treatment in the previous models was to remedy the defects of the GO approach without the scattering phase function in the calculation of scattering. However, the single-scattering albedo of the leaf is di fficult to determine in the GO approach without leaf transmittance [7]. Therefore, the single-scattering albedo of the leaf needs to adopt an empirical approach, which leads to the lack of a physical mechanism at the foundation of multiple scattering. In this study, we focused on the above limitation. In the multiple-scattering equation of row crops, we used a single-layer adding method (Figure 2) and mathematical solution of integral radiative transfer equation (its mathematical properties are also cumulative, Equation (22)) to avoid introducing the single-scattering albedo of the leaf to calculate multiple scattering. We tried to use gap probabilities to approximate the transmittance in canopy closure ( *Po*(θ*<sup>o</sup>*,*x*,*h*) and *Ps*(θ*<sup>s</sup>*,*x*,*h*) in Equation (12)) and the transmittance in the between-row area ( *Po*\_*br*(<sup>θ</sup>*<sup>o</sup>*, *h*) and *<sup>P</sup>*open in Equation (22)). In this way, the gap probabilities become a connection that can seamlessly connect single-scattering contributions and multiple-scattering contributions together.

### *4.2. Hotspot of Row Crops in the Single-Scattering Contribution*

The single-scattering contribution was accurately calculated, which was the basis of the multiple scattering calculations. Therefore, the issues of single scattering cannot be ignored in the calculation of the multiple-scattering contribution. Equations (12) and (22) require a single-scattering contribution as the original function to calculate the multiple-scattering contribution. In this study, we extended the modeling ideas of two-overlapping relationship (the overlapping relationship between leaves and canopy closures [21,22]) previously used for heterogeneity forest modeling to row modeling (Supplementary Materials B). By comparing single-scattering contributions near the hotspot in four canopies (the four growth stages of crops in Figure 4) simulated by our model and the RGM model, we found that the simulation results between the two models were very consistent (Figure 8), which indicates that our study achieved good accuracy. These results substantiate that our model considered the two-overlapping relationship to calculate gap probabilities (especially bidirectional gap probabilities and vegetation probabilities, Equations (B-24) and (B-32) in Supplementary Materials B-3), which e ffectively improved the calculation accuracy of single scattering near the hotspot in the GO approach. The mechanism for these results in Figure 8 shows that the shadow produced by the overlapping relationship between the canopy closures is gradually replaced by the shadow produced by the overlapping relationship for leaves when the row structure changes from row crops to continuous crops. In addition, the shadow area produced by the overlapping relationship between canopy closures is larger than the shadow area produced by the overlapping relationship for leaves. Therefore, the shadow near the hotspot in the crops will be reduced and single scattering near the hotspot will be enhanced. These results also imply that gap probabilities are calculated and the two-overlapping relationship is necessary.

The cause of the peak value near the hotspot involves bidirectional gap probability and bidirectional vegetation probability issues [2]. The reason for this phenomenon is that the shadow created by the overlap between the medium is not seen when the solar direction and the viewing direction are very close, and the reflectance will have a maximum value (the peak value near the hotspot) [56,57]. In continuous crops, the peak value near the hotspot is mainly influenced by LAI, leaf size, leaf orientation mode, and sun geometry [20]. However, in addition to the above factors, the peak value near the hotspot in row crops is also influenced by the row structure [31]. In our study, we found that as the row structure changes from row crops to continuous crops (from Stage\_rc1 to Stage\_cc), the width of the hotspot (slope of an inverted v-shaped area) gradually narrows, and the peak value becomes

more obvious (Figure 8 for the single-scattering contribution and Figure 6a,b for the sum of the reflectance). Our study accurately simulated the width of the hotspot, which not only considered the two-overlapping relationship but also considered the hotspot kernel function. In order to describe the width of the hotspot, cross-correlation [19] or overlap functions [57] have been commonly used as hotspot kernel functions in the RT approach. However, our study was based on the GO approach, hence the above functions were not applicable. In our study, we used the hotspot kernel function (Equation (B-25) in Supplementary Materials B-3) proposed by Chen et al. [25] in the heterogeneity forest modeling based on the GO approach. To make it suitable for the canopy of row crops, our study analogized the internal relationship between forest parameters and row parameters to the modified *C*2 in Equation (B-25), where *C*2 is the control factor of the width of the hotspot. From *C*2 (Equation (B-25)), it can be found that *L*, *h*, θ*<sup>s</sup>*, *Wp*, and *k* (note that *k* is a function of *A*1, *A*2, *h*, θ*<sup>o</sup>*, θl, and *f*(θl) in Equation (8)) were considered, which also implies that our study involved equation modeling that takes into account the main influencing factors described above, such as LAI, leaf size, leaf orientation mode, and the sun geometry, row structure. These results and equations in the modeling show that the hot kernel function is very important in simulating the peak value near the hotspot.
