## 3 Theory of Finite-Difference Time-Domain (FDTD) Method  
For the two-dimensional case studied in this work, the physical quantities in [PERSON]'s equations are independent of the z-axis. Typically, the incident electromagnetic waves are classified as either transverse magnetic (TM) or transverse electric (TE) waves. Taking TM waves as an example, according to the theory of the finite-difference time-domain (FDTD) method [29], the difference equations for TM waves in a two-dimensional electromagnetic field problem satisfy Equations (19)-(21). \[ \begin{array}{l}H_{x}^{n+1/2}(i,j+1/2)=CP(m)\cdot H_{x}^{n-1/2}(i,j+1/2)\\ -CQ(m)\cdot\frac{E^{n}(i,j+1)-E^{n}(i,j)}{\Delta y}\\ \\ H_{y}^{n+1/2}(i+1/2,j)=CP(m)\cdot H_{y}^{n-1/2}(i+1/2,j)\\ +CQ(m)\cdot\frac{E^{n}(i+1,j)-E^{n}(i,j)}{\Delta x}\\ \\ E_{x}^{n+1}(i,j)=CA(m)\cdot E_{x}^{n}(i,j)+CB(m)\cdot\\ [\frac{H_{y}^{n+1/2}(i+1/2)/\Delta x}{\Delta x}-\\ \frac{H_{x}^{n+1/2}(i,j+1/2)/\Delta x}{\Delta y}]\end{array}  \tag{19}\]  
Figure 6: The variation in the soil complex dielectric constant with moisture content. In the above equations, \(m\) denotes the position of the FDTD grid node. The parameters \(CA(m),CB(m),CP(m),CQ(m)\) satisfy Equations (22)-(25), respectively. \[CA(m)=\frac{1-\frac{\sigma(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma(m)\Delta t }{2\sigma(m)}}, \tag{22}\]  
\[CB(m)=\frac{\frac{\Delta t}{\epsilon(m)}}{1+\frac{\sigma(m)\Delta t}{2\sigma( m)}}, \tag{23}\]  
\[CP(m)=\frac{1-\frac{\sigma_{m}(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma_{m}(m) \Delta t}{2\sigma(m)}}, \tag{24}\]  
\[CQ(m)=\frac{\frac{\Delta t}{\mu(m)}}{1+\frac{\sigma_{m}(m)\Delta t}{2\sigma( m)}}. \tag{25}\]  
In introducing the incident wave at the connection boundary, a projection method is employed to determine the wave values along the boundary. Given that the projected \(y\) coordinates may fall outside the predefined sampling points within the computational domain, linear interpolation is utilized to estimate the amplitude of the incident wave at a specified distance \(d\). \[E_{inc}(d)=[1-d+floor(d)]E_{inc}[floor(d)]+[d-floor(d)]E_{inc}[floor(d)+1] \tag{26}\]  
\[H_{inc}(d)=[1-d+floor(d)]H_{inc}[floor(d)]+[d-floor(d)]H_{inc}[floor(d)+1] \tag{27}\]  
In the equation, _floor_ denotes the floor function, which rounds a value down to the nearest integer. The field value of the incident wave at the boundary is given by:  
\[E_{z}(r^{\prime})=E_{inc}(d) \tag{28}\]  
\[H_{x}(r^{\prime})=-H_{inc}(d)\cos\theta_{i} \tag{29}\]  
\[H_{y}(r^{\prime})=H_{inc}(d)\sin\theta_{i}(25) \tag{30}\]  
In the finite-difference time-domain (FDTD) region, \(\Delta x\), \(\Delta y\) represent the discrete grid spacings in \(x\), \(y\) directions, respectively. A perfectly matched layer (PML) is applied at the outer boundary of the computational domain to absorb outgoing waves. According to [PERSON]'s curl equations, the transverse magnetic (TM) mode satisfies Equations (31)-(33). \[CA(m)=\frac{1-\frac{\sigma(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma(m)\Delta t }{2\sigma(m)}}, \tag{22}\]0  
\[CA(m)=\frac{1-\frac{\sigma(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma(m)\Delta t }{2\sigma(m)}}, \tag{22}\]1  
\[CA(m)=\frac{1-\frac{\sigma(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma(m)\Delta t }{2\sigma(m)}}, \tag{22}\]2  
In the above equations, the dielectric permittivity is denoted by \(\varepsilon_{1}\), the magnetic permeability by \(\mu_{1}\), and the electrical conductivity by \(\sigma_{1}\), which are used to characterize the material properties within the computational domain. The uniaxial parameters in the \(x,y\) directions are represented by \(s_{x}\) and \(s_{y}\), respectively, and can be expressed as shown in Equation (34). \[CA(m)=\frac{1-\frac{\sigma(m)\Delta t}{2\sigma(m)}}{1+\frac{\sigma(m)\Delta t }{2\sigma(m)}}, \tag{22}\]3Among them, \(\sigma_{x}\) and \(K_{x}\) follow Equations (35) and (36), respectively.