Hybrid Method for Full-Wave Simulations of Forests at L-Band

The full-wave simulations of 3D Maxwell equations of forests at L-band are carried out with a two-step hybrid method. The tallest trees considered include one 8-meter tall trunk, 18 primary branches and 108 secondary branches. In the first step, the T-matrix of a single tree in the representation of vector cylindrical waves (VCW) is calculated by using the far-field scattering amplitudes obtained from the commercial solver FEKO. In the second step, the wave multiple scattering theory (W-MST) as represented by the Foldy-Lax equations is used to account for the wave interactions among different trees. The accuracy of the hybrid method is validated by showing that the computed solutions of three trees are in good agreement with those solved from FEKO. Next the hybrid method is used for Monte Carlo simulations of forest stands with up to 30 trees. The dependences of the solutions on tree heights and tree densities are illustrated. The results show that the full-wave based transmissivities are much larger than those obtained from the radiative transfer equation (RTE) model.


I. INTRODUCTION
Knowledge of global soil moisture distribution has a wide range of applications including agriculture management, flood assessment and climate change prediction, etc [1]. The microwave satellites such as the NASA Soil Moisture Active and Passive (SMAP) and the ESA Soil Moisture and Ocean Salinity (SMOS) have been used to monitor the global soil moisture for more than a decade and have successfully retrieved soil moisture over areas with vegetation water content (VWC) smaller than 5kg/m 2 [2]. The recent studies have reported that the data collected from the SMAP [3], [4] and the SMOS [5] are both sensitive to soil moisture over forested areas where has VWC significantly larger than 5kg/m 2 . These results encourage the studies of retrieval of soil moisture over areas with large VWC using the L-band The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . signal. A major challenge is to quantify the tree effects on the microwave signal and therefore a forward scattering model of a high fidelity is required.
The most widely used models for describing wave propagation in vegetation are the distorted Born approximation (DBA) and the radiative transfer equations (RTE). The RTE was initially proposed in [6] for modeling light transmission through foggy atmosphere and was derived heuristically based on the law of energy conservation. In the implementations for vegetation field and forests in RTE and DBA, the leaves, branches, and trunks are modeled using disks and cylinders of difference sizes. A significant number of adhoc assumptions are made in the RTE which are listed as follow: (i) The positions of branches and leaves are assumed to be statistically homogeneous in space. (ii) The vegetation canopy is modeled as an effective medium with an effective propagation constant and an associated effective opacity ''tau''. (iii) The scatterers such as trunks, branches and leaves are in the far field of each other and they are also in the far field from the soil below. Several analytical efforts have been made to improve the accuracy of RTE during the last few decades [7], [8], [9], [10], [11], [12], [13], [14], [15]. The coherent additions of scattering amplitudes from branches and leaves was applied to investigate the effects of the plant structures [12]. The row structures within the corn field are studied using the periodic slap together with the phase matrix [11]. The coherent model is based on superposition of first order scattered fields. It is well known that the first scattering violates energy conservation. Second order scattering has to be included so that the combined first order and second order obey energy conservation. The theorem can be shown analytically by perturbation theory. For explicit expressions they have been shown in [16] for rough surface scattering and [17] for point scatterers in volume scattering.
With the advancement in the digital computers, full-wave simulations of the vegetation have been attempted using the numerical method such as the method of moments (MOM) [18], [19], the finite element method (FEM) and the finite-different time-domain (FDTD) [20]. However, due to the tremendous memory and computation cost, the simulations are limited to several plants which are unable to represent the real vegetation field. Since 2017, we have been developing a two-step hybrid method for performing full-wave simulations of a large vegetation field [21], [22], [23], [24]. Our method consists of (i) treating each plant or each tree as one scatterer instead of a single branch or a single leaf as scatterer, (ii) each plant or tree is enclosed by a circumscribing vertical cylinder, and (iii) vector cylindrical wave expansions (VCW) are used for the scattered fields outside the circumscribing cylinders. In the first step of the hybrid method, the T-matrix of a single plant or a single tree is calculated. In the past, T-matrices were based on vector spherical waves (VSW) and the extended boundary conditions (EBC) were used to calculate the T-matrix. This method is not appropriate for plants as vector cylindrical waves (VCW) are more suitable for representing them than VSW. The EBC method, using complete domain basis functions, is unstable for objects of large aspect ratio. In our hybrid method, the commercial softwares such as the MOM based FEKO and the FEM based HFSS are utilized to perform full-wave simulations of one single plant, from which the corresponding T-matrix is extracted based on the near-field using the Huygens principle and VCW expansions [23]. The full-wave based T-matrix characterizes the scattering of one single plant and captures the multiple scattering effects caused by the plant structures. In the second step, the T-matrix is combined with the Foldy-Lax multiple scattering equations (FLE) to consider the wave interactions among different plants. The resulting closed-form equations are derived analytically using the VCW expansions and the translational addition theorem. The full-wave simulations of a large wheat field that has up to 169 plants at the L-, S-, and C-bands were performed using the hybrid method to investigate the vegetation effects on the microwave [24]. The results show that the full-wave based microwave transmissivity is much larger that of the RTE and the absorption cased by the vegetation will saturate as frequency increases.
Driven by the needs of mapping soil moisture over forested area, we applied the hybrid method to simulate the forest in this paper. The problems of simulating one tree as well as a forest stand are much more challenging than simulating other plants and crops because trees are much taller. In our previous paper [22], we considered the forest canopy modelled by 196 vertical dielectric cylinders in full wave simulations at L-band. Because only the trunks were considered, the circumscribing cylinder had a radius of only R c = 6cm. In this paper, we add a significant number of primary branches and secondary branches which extend out from the tree trunk. Thus, the circumscribing cylinder of a much larger radius, such as R c = 1m which is 16.6 times larger than that of the previous paper, is needed to enclose the entire tree. The increase in size creates two challenges for the hybrid method. The first challenge is the matrix inversion-based approach [23] for calculating the T-matrix of single plant fails when the plant has a large electrical size and a complicated structure. This method requires a careful selection of the incident angles such that all relevant information of the scatterer are captured within the corresponding scattering coefficients. In addition, the number of incident waves to be simulated is equal to the number of the modes used in the T-matrix. The reason that the matrix inversion method works for wheat is because a simple stem only structure is used to model the wheat and thus a uniform sampling of the incident angle is good enough for the T-matrix extraction. To overcome this issue, a general relation is discovered in this paper which shows that T-matrix is the coefficient of the Fourier series of the scattered field coefficient. Based on that, a far-field based approach is developed for calculating the T-matrix of one single tree.
The second challenge in simulating a forest stand is the memory demand for solving the closed-form FLE as it increases dramatically with the number of modes used in the T-matrix. This issue is resolved by adopting a physical iterative solution for considering the multiple scattering among different trees, which eliminates the need for the direct solution of FLE. First, the accuracies of the two new proposed improvements in the methodology are validated by comparing the results with FEKO for solving scattering from three 8-meter tall trees. Then, the full-wave Monte Carlo simulations of forests with different heights and densities are performed using the hybrid method to investigate their impacts on the microwave transmissivity. The results are compared with those of the RTE model. In this paper, the tallest trees considered are 8 meters tall and include the trunk, 18 primary branches and 108 secondary branches.
The paper is organized as follow. In Section II, we discuss the hybrid method, which includes the forest configuration and the tree structure, the derivation of the far-field based T-matrix, the procedure of physical iterative solution for solving FLE, and the plane wave transformation of VCW for inner field calculations. In Section III, the hybrid method is validated with FEKO by solving scattering from three 8-meter tall trees. In Section IV, the hybrid method is applied to perform full-wave simulations of forests and the resulting microwave transmissivities are compared with the RTE. The conclusions are drawn in Section V.

II. HYBRID METHOD A. FOREST CONFIGURATION AND TREE STRUCTURE
The hybrid method is built based on a realistic forest stand configuration that can capture both structures of trees and the gaps among different trees. The top view of the forest to be simulated is shown in Fig. 1, which is a circular area of radius R. Each small circle denotes one single tree and its radius R c is the smallest circle that can fully enclose the corresponding tree. For the tree considered in this article, the R c equals to the length of the primary branch as listed in Table 1. The positions of the trees are randomly generated within the circular area and the number of trees is determined by a density parameter. Since each tree is treated as one single scatterer in the hybrid method, we should make sure that they do not overlap with each other when generating their positions. The components of one single tree include trunk, branches of different sizes and leaves. The full-wave simulation of the entire tree is still a challenge as the number of leaves can be in a million-scale and thus cannot be done given the available computation resources. In this paper, we do not consider the small branches and the leaves in the tree model. Specifically, a tree is represented using the trunk, the primary and the secondary branches as shown in Fig. 2. Although a realistic structure can be used in the hybrid method as a numerical solver is employed to perform the single plant fullwave simulations, all the components are modelled using the dielectric cylinders of different sizes to make the result comparable to the classical RTE model. In the future, the realistic tree structure can be used to further improve the simulation. The full-wave simulations of the branch only tree can be done using the Multilevel Fast Multipole (MLFMM) solver of FEKO. Based on the forest stand configuration, the full-wave simulations of the stand are performed using the hybrid method in two steps.

B. T-MATRIX EXTRACTION OF ONE SINGLE TREE 1) DEFINITION OF T-MATRIX
The T-matrix is a general approach for describing the response of one single scatterer and is first proposed in [25] based on the spherical wave. The T-matrix concept is extended to cylindrical coordinate in our research to characterize the response of one single plant based on the VCW expansions using the full-wave solutions obtained from the numerical solver. To construct the T-matrix, an infinite cylindrical surface is used to enclose a single plant such that the space is divided into two regions. While the plant is in the inner region, the outer region is the empty space as shown in Fig. 2. Based on the model expansion theory, all the field in the cylindrical coordinate can be expanded in terms of VCW. Specifically, the scattered field in the outer region can be expressed as whereM n (k z ,r) andN n (k z ,r) are the outgoing VCW, a SM n (k z ) and a SN n (k z ) are their corresponding expansion coefficients. The outgoing VCW are defined as [16] wherer =ρρ +φφ +ẑz is the position vector in cylindrical coordinate, k ρ = k 2 − k 2 z is the propagation constant in theρ direction, n is the order of harmonic inφ direction, the H (1) n (k ρ ρ) is the Hankel function of first kind and H (1) n (k ρ ρ) is the corresponding first order derivative. The excitation field in both regions can be expanded as where RgM n (k z ,r) and RgN n (k z ,r) are the incoming VCW and their expressions are obtained, respectively, by replacing the H (2) and (3) with Bessel function J n (k ρ ρ). The corresponding excitation coefficients are denoted respectively as a EM n (k z ) and a EN n (k z ). They are related to the scattered field coefficients through the T-matrix This is the definition of the VCW based T-matrix. If the T-matrix of one scatterer is known and the excitation is given, then the corresponding scattered field coefficients can be calculated from (5) and (6). Then the resulting scattered field is obtained using (1). An advantage of the hybrid method is that the T-matrices calculated in step 1 are re-usable. Using a calculated T-matrix, the T-matrix can be rotated to give the T-matrix of a different tree because a tree does not have azimuthal symmetry. By calculating several trees with varying heights, a library of T-matrices can be generated. The library of T-matrices can be used for many realizations of varying positions and varying rotations in a forest. Thus the hybrid method also works for forest stand with mixed trees.

2) T-MATRIX AND SCATTERED FIELD COEFFICIENT
Let the excitation be a plane wave given as where k is the propagation constant,k i is the propagation vector,v i andĥ i , respectively, are the vertical and horizontal polarization vector of the incident direction (θ i , φ i ), E vi and E hi are their corresponding amplitudes. The triplet (k i ,v i ,ĥ i ) forms an orthonormal unit system and are defined aŝ The VCW coefficients of the plane wave are known as [16] a EM n (k z ) = E hi where δ is the Dirac delta function. Substitutions of (9) into (5) and (6) and the cancelation of integration over k z with the delta function give Note that the two cylindrical wavesM n (k z ,r) andN n (k z ,r) correspond to the H-pol and V-pol respectively. By applying an H-pol (E vi = 0, E hi = 1) and a V-pol (E vi = 1, E hi = 0) plane wave, respectively, to (10) and (11), the Fourier series of the scattered field coefficients can be obtained where β = M , N denotes the polarization, the prime superscript indicates the source, and . It reveals that the T-matrix element of T ββ nn (k z , k zi ) is the coefficient of the Fourier series of the plane wave scattered field coefficient α Sβ n,β (k z , k zi ) in terms of φ i . Multiplying both sides with dφ i exp(in φ i ) and using the orthogonality of the exponential function the T-matrix can be calculated using the scattered field coefficient from The significance of (14) is that it enables us to compute the T-matrix element by element and removes the inherent constrain set by the matrix inversion-based approach for calculating the T-matrix [26].

3) FAR-FIELD BASED T-MATRIX EXTRACTION
It is important to notice that the expansion coefficients in (1) are independent of distance. Given the a Sβ n (k z ), the corresponding scattered field are known regardless of the distance. Conversely, a Sβ n (k z ) can be obtained either from the near-field or the far-field. The near-field based approach for calculating the a Sβ n (k z ) using the Huygens principle and the VCW expansions of the Dyadic Green's function has been presented in [23]. In this paper, the far-field based method for extracting the T-matrix is developed using (14).
In the far-field region, the Hankel function can be approximated as Substitute (15) into (2) and (3), we can get the far-field approximations of the VCW lim ρ→∞M n (k z ,r) (17) where theM n (k z ,r) only hasφ component and theN n (k z ,r) hasρ andẑ components. With the far-field approximations, the scattered field in the far-field region becomē By applying the stationary phase approximation to (18), the far-field can be simplified as Note thatv =ρ cos θ −ẑ sin θ andφ =ĥ, the scattered farfield becomē E s (r) = −k sin θ n va SN n (k cos θ) + iĥa SM n (k cos θ) It is also noticed that the last term is the far-field approximation of the free space Green's function which is known as [16] e ikr r = Using (21), the far-field is related to the VCW coefficient through In the far-field, the scattered field is related to the incident field through the scattering amplitude f ββ (k s , k i ) By equating (22) and (23), the scattered field coefficients are related to the far-field amplitude 2k sin θ 2k sin θ n i 1−n exp(inφ)a SM n (k cos θ)(i) Applying an H-pol (E vi = 0, E hi = 1) and a V-pol (E vi = 1, E hi = 0) plane wave excitation respectively to (24) and (25), and integrating both sides with exp(inφ), we can obtain a Sβ n,β (k cos θ) = where . This states that the scattered field coefficient a Sβ n,β (k cos θ) is the coefficient of Fourier series of the farfield amplitude f ββ (k s , k i ) in term of φ. Combination of (14) and (26) gives (27) It indicates that the T-matrix elements are the Fourier coefficients of the double Fourier series of the far-field in terms of φ and φ i . Consequently, we can use the far-field to extract the T-matrix.

C. ITERATIVE SOLUTION OF FL EQUATIONS
The second step of the hybrid method is to combine the T-matrix with the FLE to account for the multiple scattering among different trees through where N t is the number of trees, G qp is the free space Green's function representing the propagation from tree p to q, and theĒ ex q (r) denotes the final excitation acting on the q-th tree. The FLE states that the final excitation fieldĒ ex q (r) acting on the q-th tree is the incident field plus the final scattered fieldĒ sc p (r) of all other trees except itself. In (28), the final scattered field is expressed in terms of T-matrix as E sc p (r) = T pĒ ex p (r) and thus a set of self-consistent equations are formed withĒ ex p (r) being the only unknowns. The FLE are derived from the Maxwell equations [27]. The numerical solutions of FLE are the exact solutions of Maxwell equations. The exact fields are governed self-consistently by the Maxwell equations and boundary conditions on the scatterers.
By using the VCW expansions and the translational addition theorem, a set of linear equations are obtained for solving the final unknown excitation coefficient. Previously, the closed-form equations were solved directly using the numerical method of GMRES to obtain the final solution [24]. However, this becomes impractical for forest simulation as the number of unknowns is on the order of O(N 2 t N 2 T ), where N T is the number of modes used in the T-matrix. To avoid the memory challenge, the FLE can be solved iteratively based on the physical scattering process such that we do not need to store the FLE matrix. The iterative process starts from the first order scattered field which is excited by the incident field only and can be obtained from E s (1) q (r) = T qĒ inc (r) (29) where the order in the superscript indicates the number of time that the field is scattered. The higher order scattered field is then obtained from the lower order using the recurrent iterationĒ which states that the total exciting field acting on q-th scatterer is the sum of the (n − 1)-th scattered field from all other scatterers. The scattered field will decay as the order increases and the iteration can be terminated when the higher order scattered field is negligible. The total scattered field is then calculated fromĒ where N it is the total number of the physical iterations. The iterative solution is preferable when the density of scatterers is small and the solution will converge in a few iterations as we will show in the result section.

D. INNER FIELD EXTRACTION
The solutions solved from the FLE are the final excitation field coefficients of each plant a Eβ n,q (k z ). The corresponding final scattered field coefficient a Sβ n,q (k z ) can be obtained from T-matrix through (5) and (6). If the observation point is in the outer region, thenĒ s q (r) can be computed from (1) using a Sβ n,q (k z ). It is important to notice that the T-matrix is only valid forĒ s q (r) in the outer region and thus can not be applied to compute the inner region scattered field. The scattered field of each plant in its inner region can be obtained using the plane wave transformations and the superposition principle. The solutions solved from the FLE express the final excitation field acting on each scatterer in terms of VCW using (4). By applying the plane wave expansions of the VCW which are known as [16] RgM n (k z ,r) we can transform (4) intō +C v q (θ k , φ k )v (θ k , φ k ) exp(inφ k ) exp(ik ·r) (34) which expresses the final excitation field as a linear combination of the plane waves from different angles (θ k , φ k ), C h q and C v q are the corresponding plane wave coefficients and they are related to the VCW coefficients as Using the linearity of the electromagnetic field and the superposition principle, the corresponding inner region scattered field excited by (34) can be obtained from whereĒ s,β q,in (r) is the scattered field of the tree excited by the β-pol plane wave from (θ k , φ k ) and has been already computed in the first step when calculating the T-matrix. After computing the scattered field from each tree, the scattered field from the whole forest stand is achieved bȳ E s (r) = N t q=1Ē s q (r). Then the total field on the bottom plane can be calculated from which can be used to compute the microwave transmissivity.

III. METHOD VALIDATION
In this section, we validate the far-field based approach by calculating the T-matrix of a 8-meter tall tree and the accuracy of the hybrid method is tested with FEKO by solving the scattering from three 8-meter tall trees. The excitation field is a vertically (V-pol) or horizontally (H-pol) polarized plane wave from the direction of φ i = 0 o , θ i = 40 o , at L-band 1.41 GHz, which are the incident angle and the operation frequency of NASA Soil Moisture Active Passive (SMAP) mission [28]. Three types of components, the trunk and, the primary and secondary branches, are considered in the tree model. Their sizes and numbers are listed in Table 1.
While the trunk is a vertical cylinder, the orientation angles (φ B 1 , θ B 1 ) of the primary branches are randomly generated with a uniform distribution of . Three 8-meter tall trees generated based on these parameters are shown in Fig. 4. A constant relative permittivity r = 13.68 + i4.62 obtained from the dielectric model presented in [29] is used in the simulation. The corresponding volume fraction of water is m v = 20%.

A. T-MATRIX
First, we perform full-wave simulations of an 8-meter tall tree using the FEKO to obtain the scattered far-field from which the corresponding T-matrix is calculated using (27). The VCW based T-matrix has two free variables n and k z which range from −∞ to +∞ in theory. In practice, they are truncated at a maximum value that ensures the convergence of the scattered field. The maximum number of modes in theφ direction is related to the R c as N max n = 2kR c . The continuous integration over k z is truncated at k max z and is approximated using a N θ uniform sampling within [−k max z , k max z ]. In this paper, the T-matrix only consider propagating waves as the evanescent waves are negligible so k max z = k 0 . The number of samplings N θ is chose based on the convergence test. Consequently, the total number of modes used in the T-matrix is N T = 2(2kR c + 1)N θ . The far-field based T-matrix (27) requires a continuous integration over φ i which is very expensive as it needs the full-wave simulation for each sampled φ i . To improve the simulation efficiency, the Legendre quadrature is applied to the φ i integration and the resulting number of samplings is N φ i = 2n which is 5 time less than 10n which is the sampling number required by a uniform sampling.
To justify the truncation we made in the number of modes and validate the accuracy of the T-matrix, the scattered field of an 8-meter tall tree excited by a V-pol plane wave on a vertical z-line is plotted in Fig. 3(a). The amplitude of incoming wave is a constant |Ē inc | = 1V /m. We can see despite a small mismatch on one of the peak caused by the sharp scattering from branch, the scattered field obtained from the T-matrix agree very closely with the FEKO solution. As the evanescent waves are largely caused by the tree trunk and decay exponentially in theρ direction, they are negligible on the enclosing cylindrical surface where ρ = 1m and is about 4λ away from the trunk. That is the reason why we can exclude the evanescent waves in the T-matrix. The scattered field over a horizontal x-line is also plotted in Fig. 3(b) which shows that the scattered field decay as ρ increases. This implies that when the distance between two trees is greater than a certain value ρ max , their interactions are negligible. Consequently, the area of the forest stand to be simulated can be truncated at R = ρ max . That is the reason for the simulated forest stand is a circular area of radius R = 8m. The good agreements of the field solutions validate the accuracy of the far-field based T-matrix.

B. SCATTERING FROM 3 TREES WITH VALIDATION USING FEKO
Second, we validate the iterative solutions of the FLE by solving the scattering from three 8-meter tall trees as shown in Fig. 4. The iterative solutions of the scattered E-field along the z-axis from −4m to 4m are plotted in Fig. 5. The amplitude of the scattered field decreases dramatically with the order of scattering and is on the order of 10 −6 at the 10-th iteration. Thus the iterative process can be terminated at N it = 10. The hybrid method solutions are obtained by summing up the 10 order of solutions and are compared with the FKEO solutions in Fig. 6. The good agreements of the near-field solutions validate the accuracy of the physical iterative method. Since the hybrid method treats each single tree composed of many components, rather than a single branch or trunk, as an individual scatterer and the typical density of tree ρ t in the forest is within [0.07, 0.15](1/m 2 ), the number of scatterers (tree) is always small. For example, the number of trees within a circular area of R = 8m with ρ t = 0.15 is 30, which is much smaller than that of the crop. Thus, a fast convergence can always be achieved by the physical iterative process when simulating the forest. Fig. 7 shows the scattered E-field on a plane 5 cm below the trees (see Fig. 4) calculated from the two methods. Despite a small discontinuity between the fields in the inner and the outer regions, the hybrid method results match well with those calculated from the FEKO. This validates the wave transformation approach for computing the inner region field. Moreover, the coherent wave approach enables the hybrid method to model the field variations caused by the multiple scattering due to both the tree structures and spatial variations of the tree positions. Consequently, the non-uniform transmission within the forest can be captured and is much  more accurate than the classical incoherent RTE which only uses a constant extinction coefficient κ e to characterize the attenuation caused by the vegetation layer.

IV. FULL-WAVE SIMULATIONS OF FORESTS
In this section, we apply the hybrid method to calculate the microwave transmissivity at L-band over the forested area. We perform full-wave numerical solutions of Maxwell equations in Monte Carlo simulations of trees in a physical area of 200 m 2 . We vary tree densities and tree heights to investigate their impacts on the microwave transmissivity. The top view of the forest is shown in Fig. 1 with a radius R = 8m. The T-matrix of a single tree is first calculated as described in Section II based on the full-wave solutions obtained from the FKEO. The T-matrix is stored and can be repeatedly used in different realizations of the Monte Carlo simulations. The reusability of the T-matrices is an important advantage of the hybrid method. For each realization, the positions of the VOLUME 10, 2022 trees are randomly generated within the circular area to form a different realization of forest. Based on that, the hybrid method is applied to perform the full-wave simulations and the total field of the forest stand can be obtained following the procedure described in Section II-D. Then the transmissivity of the microwave is obtained by calculating the normalized transmitted power over a reception area through whereS is the Poynting vector and A is a square reception area of size l x at the bottom of the forest as shown in Fig. 1. Note that the reception area is different from the physical area and is smaller than the physical area to avoid edge the effects caused by the truncation of forest. The final transmissivity is obtained by the ensemble average over T computed from N r realizations using In the following, we first perform the convergence tests to determine the numerical parameters that are used in the Monte Carlo simulations. The convergence tests are unique for the random media scattering problem and they are not required in the deterministic problem. They test the convergence of t with respect to the reception area A and the number of realizations N r . After that, we will apply the hybrid method to study the effects of the tree heights and the tree densities on the microwave transmission in the forested area.

A. CONVERGENCE TESTS
Wave propagation in the forest is a random scattering process and the microwave transmissivity T is a random variable.
In microwave remote sensing application, only the mean value of T is of interest and is computed using (40) based on the Monte Carlo simulations. Generally, the randomness of T is proportional to the tree height, thus the convergence tests are performed over the tallest forest (H = 8m) considered in this paper to determine the numerical parameters of l x and N r . We first test convergence of transmissivity over the size of the reception area A. The T as defined in (39) is a normalized value that will converge to a value for a large enough area A. Since a limited size of forest is simulated, the reception area is thus located at the center region as shown in Fig. 1 to avoid the edge effects caused by the truncation of the forest. Fig. 8 shows the L-band H-pol transmissivities of the 8-meter tall forest stand as the function of A. The results are calculated based on 20 realizations. We can see the t of different tree densities converge as the area increases indicating that the simulated area is large enough to obtain a converged transmissivity. Considering that we need to minimize the edge effects, it is preferable to keep A as small as possible while maintaining the convergence of t. Thus, the size of the reception area is chosen to be A = 4m 2 with l x = 2m.
Next, we carry out the convergence tests of t over the number of realizations N r . Fig. 9 shows the transmissivities of three different tree densities are plotted against N r . While the t of the two sparser cases converge at N r = 10, it takes 7 more realizations for the denser case to converge. This makes sense as a denser forest has stronger multiple scattering, thus results a larger fluctuation over T obtained from different realizations. Consequently, it requires more realizations for a denser forest to obtain the converged t. Based on the testing results, the number of realizations used in the following simulations is N r = 20.

B. TRANSMISSIVITY VS. TREE HEIGHTS
After determining the numerical parameters, we perform Monte Carlo simulations of forests with three different heights to investigate its impact on the microwave transmissitivty. The positions of the trees are generated based on the 8-meter tall trees and the same positions are used for the forest stands of three different heights as listed in the Table 1.
To demonstrate the importance of full-wave simulations, the full-wave calculated transmissivities are compared with those obtained from the RTE using the same parameters. In RTE,  the scatterers are the individual trunks, primary branches and secondary branches. The positions of these scatterers are assumed to be statistically homogeneous within the forest. Thus, the forest canopy is approximated as an effective layered medium with an effective attenuation coefficient κ e and the transmissivity is estimated from where d is the thickness of the vegetation layer and equal to the tree height, θ i is the incident angle. The κ e describes the power attenuation rate per unit distance in the forest and is computed from where n 0 j is number of j-th component per m 3 , σ s j and σ a i are the corresponding scattering and absorption cross sections, and · is the ensemble average on the orientation angle (θ c j , φ c j ). Since all the components are modelled as dielectric cylinders, their σ s j and σ a j are calculated analytically based on the ICA using the optical theorem [16]. The transmissivities of different tree heights calculated from the two methods including both the V-pol and H-pol results are plotted in Fig. 10. While a similar trend that t decreases with H is observed for both methods, the hybrid method results are significantly larger than those of the RTE. For instance, the full-wave based V-pol transmissivity for 8-meter tall forest is about 0.59 which is 55% higher than 0.38 predicted by the RTE. As the approximations made in the RTE become more inaccurate as the tree height increases, their difference also becomes larger as tree becomes higher.
C. TRANSMISSIVITY VS. TREE DENSITIES Next, we simulate forests of different tree densities to investigate the dependence of transmissivity on tree density. The simulations are performed for forest stands of 5m-and 8m-tall trees using both the V-pol and H-pol plane wave. The tree density varies from 0.07 to 0.15 and the corresponding number of tree ranges from 14 to 30. The microwave transmissivities obtained from the hybrid method and the RTE model are plotted in Fig.11. The two methods both predict that the H-pol transmission is greater that of the V-pol. This is physically due to the fact that the attenuation is predominantly caused by the vertical trunk. The full-wave based transmissivities are significantly greater than those of the RTE. This is because the homogenization approximation made in the RTE is only valid for a dense vegetation canopy and is not applicable for forest whose scatterers has a clustered spatial variation. In RTE, an effective κ e is introduced in an adhoc manner to characterize the wave propagation which do not follow from Maxwell equations. Such adhoc assumption also implies that there is constant attenuation in all the space of the vegetation layer. However, this contradicts to the fact that there is no attenuation in the free space and only the vegetation cause the absorption. The failure of RTE in capturing the gap spacing results in an underestimation of the microwave transmission. The full-wave approach can model the free space propagation as it includes both the amplitude and phase of the scattered waves from each tree.

V. CONCLUSION
In this paper, we presented Numerical solutions of Maxwell equations in 3D (NMM3D) of forests at L-band with trees that are up to 8-meter tall. A realistic tree model, including the tree trunk, the primary and secondary branches, were used. Simulations were performed for tree densities up to 0.15 per square meter and a maximum of 30 trees were used. Results show that full-wave simulations give significantly higher transmissivity than RTE. Improvements of methodology that are made in this paper include the use of the scattering amplitudes obtained from FEKO to calculate T-matrix coefficients and the use of physical iterative approaches for solving the FLE. In the future, the leaves and twigs will be included to better represent the tree and a new approach will be developed for calculating the T-matrix. In 2022, the SMAPVEX campaign have been carried out over the forested area to measure both the ground true and the microwave responses at L-and S-bands. The measured data can be used to further validate the results. Presently we are studying the inclusion of evanescent waves for the T-matrix in the VCW representation. The numerical solutions of the hybrid method will be accelerated. A fast method for solving the FLE is also being developed and a comprehensive comparison of computation cost between the fast hybrid method and the classical numerical solver will be made in a future publication. With significant improvement in computational efficiency, comparisons with experiments will be made.