Time-Domain Electromagnetic Wave Propagation in Confined Environments

Confined environments like tunnels are electrically large structures for guided wave propagation. They can have arbitrary cross sections, and the design and optimization of antenna for communication system requires the knowledge of a “full-wave” solution in nearby zones. Current models based on asymptotic approaches do not describe adequately the wave propagation under the above conditions. In addition, a complete “full-wave” analysis of the tunnel propagation performances is not feasible in terms of computer expenditure. After a survey of the most commonly used techniques for propagation in tunnels, some investigation regarding an appropriate approach to find the fields is proposed. It is based on a modal decomposition of the wave propagation that allows an optimization of the coupling with the antenna. To find the mode characteristic for arbitrary cross section, a full-wave method, namely, the transmission-line matrix (TLM), is modified to a so-called 2.5-dimensional TLM algorithm and presented in details. This approach is validated for a canonical structure. Then, it is applied to study the wave propagation in a realistic rectangular tunnel. The concept of surface impedance boundary condition (SIBC) is introduced to reduce the TLM computational domain and model the tunnel walls that can be considered as lossy dielectric. Results show that guided structures with lossy dielectric walls of arbitrary cross section can be studied with this approach.


Introduction
Railway, roadway and mine tunnels, buildings, and warehouses are some examples of confined environments, in which electromagnetic wave propagation has to be investigated With the increasing development of computers, appropriate new models and simplifications are being developed. The formulation of an efficient modal approach stemming from the fact that tunnels can be modeled by an over-sized waveguide and the a priori knowledge of the fields in the axial direction is presented. It allows one to have a better physical insight into wave propagation in confined environments, as well as dealing with electrically large structures like tunnels. The mode parameter determination is carried out by a full-wave time-domain method, namely, the transmission-line matrix (TLM) method.
The calculation volume has to be limited in full-wave volumic methods, and we are interested only in the fields inside the tunnel. Thus, the electromagnetic modeling of the tunnel walls, which can be lossy dielectric, is addressed. Lastly, by assimilating a confined environment to a lossy dielectric waveguides of arbitrary cross section, the mode extraction of these structures is presented. To the best of the authors' knowledge, no such model has been reported.
The chapter is structured as follows: In the first section, an overview of the principal techniques for the description of the EM wave propagation in above structures is briefly presented. In the second section, the formulation of the modal approach is described in detail. In the following section, the implementation for a simple canonical case is shown. Lastly, the numerical analysis of multimodal waveguides representative of confined environments is illustrated in a realistic rectangular tunnel. Finally, discussions and conclusions are developed at the end of the chapter.

Modeling approaches for wave propagation in confined environments
The correct understanding of wave propagation in confined structures like tunnels has been an important area of research and development. They have been studied in the last 40 years for radio communication system deployment. Unfortunately, as we shall see, current models cannot be generally applied or they do not describe adequately the wave propagation. The ❚✐♠❡✲❉♦♠❛✐♥ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡ Pr♦♣❛❣❛t✐♦♥ ✐♥ ❈♦♥✜♥❡❞ ❊♥✈✐r♦♥♠❡♥ts ✸ ✶✵✳✺✼✼✷✴✻✶✸✹✵ purpose of this section is to provide an overview of the approaches for the problem of wave propagation in these structures.

Analytical models: circular and elliptical dielectric tunnels
No analytic solution exists for the problem of wave propagation along a dielectric waveguide of arbitrary cross-sectional shape. The circular and elliptical waveguides are the only cases for which analytic solutions for guided waves exist. Field solution in waveguides for analytical and approximate methods is usually expressed in terms of modes. They are classified into different types according to their field configuration: transverse electric modes (TE), with no electric field in the direction of propagation; transverse magnetic modes (TM), with no magnetic field in this direction; transverse electromagnetic modes (TEM), with neither electric nor magnetic field in this direction; and hybrid modes (HE) or (EH), with nonzero electric and magnetic fields in the direction of propagation.
Dielectric circular waveguides can support a family of circularly symmetric TE nm or TM nm modes, with n = m = 0, and hybrid modes HE nm and EH nm . The subscripts n and m denote the number of oscillations with the cylindrical coordinates ρ and φ, respectively. Some discussions considering different excitations, dependence on the constitutive parameters of the tunnel walls, and other interesting topics are treated in [9]. The details of the determination of the fields will not be repeated here. Finally, the wave propagation in elliptical tunnels has not been treated. However, a deformed circular waveguide can be approximated by an elliptical cylinder. Detailed theoretical as well as experimental results can be found in [34].

Approximate models: rectangular dielectric tunnels
An exact analytic solution does not exist for the case of wave propagation in a rectangular tunnel with lossy dielectric walls. An approximate approach, namely, Marcatilli's method, is usually employed to analyze this structure. The a priori assumption of Marcatilli's approach is that most of the energy of the modes propagating in the structure is contained within the core region and very little guided power is contained in the corner regions of the guide [34]. Then, the boundary conditions are only matched along the four sides of the hollow region, and the fields in the corners are ignored. The Helmholtz equation in Cartesian coordinates is solved with separated variables, and approximate solutions for the propagation constant and for the field distributions are found. The modes are classified into E x nm and E y nm with most of its electric field polarized in the horizontal or vertical direction, respectively. This approximation gives the mode parameters with sufficient accuracy as long as the assumptions below are valid: ✹ ❆❞✈❛♥❝❡❞ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡s where n and m are integers describing the nm-th mode, λ is the free-space wavelength, w and h the waveguide dimensions, andǭ is the normalized complex permittivityǭ = ǫ r − jσ/ (ωǫ 0 ), where σ is the wall conductivity. Finally, a study of modal propagation in curved rectangular tunnels can be found in [20].

Asymptotic models
Ray models use high-frequency expansions of Maxwell's equations and are based on optical laws to model the fields. They constitute the most commonly used deterministic techniques to analyze confined environments [32]. In ray-based methods, the waves are treated as rays propagating perpendicular to the wave fronts. This approximation is valid at positions sufficiently distant from the source. The resulting spherical waves can be approximated by a locally plane wave on a small portion of the sphere. Two methods are widely used: ray launching and ray tracing [8,10]. Curved tunnels, arched cross sections, or bends of a tunnel can be treated by ray methods provided that radii of curvature of the surface be large compared to the wavelength [22,33].
Current tunnel cross sections are getting smaller, and antennas have to operate close to the tunnel walls, such that these interfere with their near-fields. As a result, antenna performances cannot be guaranteed by using asymptotic methods. For this reason, full-wave models accounting for near-field effects become more relevant.

Full-wave models
In electromagnetism, two main families of numerical methods exist: frequency and time-domain methods. Nowadays, time-domain methods are becoming increasingly used [11,31]. They are attractive because of their relative simplicity and their ability to account for arbitrary geometries, obstacles like vehicles, nonlinearities, and different materials and to determine transient response and perform wideband characterization [6,32,35]. Thus, full-wave methods attract more attention for the purpose of deterministic radio coverage predictions, near-field considerations, simplifications, and practical implementations [13]. However, efficient and adequate models are still needed to describe the wave propagation in tunnel environments. A modal approach is preferred because it allows one to describe the propagation by modes [1,2,29]. An efficient time-domain modal approach based on the transmission-line matrix (TLM) method will be introduced in the next section.

The Transmission-Line Matrix (TLM) method for guiding structures
The transmission-line matrix (TLM) method was developed by P. B. Johns and his co-workers in the 1970s [17]. Its theoretical foundations are based on Huygens's model of wave propagation. The TLM method is very attractive and flexible to analyze electromagnetic field problems and has received increased recognition for full-wave analysis of arbitrary-shaped guiding structures [15]. Through the years, several nodes have been developed in TLM for structured and unstructured meshes in two and three dimensions. We will focus on the symmetrical condensed node (SCN), which is the most widely used for 3D structures.
Advanced Electromagnetic Waves 242

Description of 3D TLM method
The symmetrical condensed node (SCN) consists of a network of interconnected multi-port devices where the electromagnetic wave propagation is simulated by the propagation of traveling pulses; a unit cell is presented in (Figure 1). On each arm, two orthogonal ports are employed to account for any polarization. To solve an electromagnetic problem, the solution region has to be divided into a number of these elementary nodes.
A rule of thumb consists of taking the maximum cubic cell size as ∆l = λ/10, where λ correspond to the medium wavelength of the highest frequency of interest. Then, excitation is imposed. The objective of the source terms is to simulate a desired phenomenon in the structure. For instance, in guiding structures, a given mode can be selectively excited by an adequate field distribution that corresponds to the mode transverse configuration (mode template). Some other problems require to analyze the structure over a wide frequency range. This is achieved by exciting the structure with a wideband time signal such as Dirac, Gaussian, etc. Thus, in general, the excitation corresponds to a space-time distribution. Nodes are interconnected by these virtual transmission lines, and the excitation propagates from the source nodes to the adjacent nodes at each time step.
The method is carried out basically through two processes: scattering and connection. In the scattering process, voltage pulses n V i are incident upon the node from each of the link-lines (halfway between two nodes) at each time step n∆t. These pulses are then scattered to produce a set of scattered voltages, k V r , which become incident on adjacent nodes at the next time step (n + 1)∆t. In the connection process, pulses are simply exchanged among immediate neighbors.
Volumic methods such as TLM full-wave methods require a limitation of the computational domain for open problems in which fields exist at large distances. Thus, somehow the problem has to be bounded. The boundary conditions are the set of conditions which specifies the behavior of fields at the boundaries of the computational domain.

Boundary conditions in TLM
The boundary conditions link the electromagnetic fields through the tangential or normal field values. Since TLM is based on the equivalence between Maxwell's equations and equations for voltages and currents that travel in a mesh of interconnected transmission lines, a relationship of the involved voltages at the boundary can be found. The scattered voltages V r are always known values, and the incident voltages V i are unknown. Any resistive load at a boundary may be simulated by introducing a reflection coefficient Γ as shown in (3), [7].
This formalism allows us to represent a variety of boundary conditions as long as a reflection coefficient Γ can be defined. For instance, for a perfectly electric conductor (PEC), boundary is simulated by choosing Γ = −1; a perfectly magnetic conductor (PMC) is implemented by choosing Γ = 1. The reflection coefficient for lossy boundaries [15], relating the incident and reflected voltages, can be expressed in Laplace domain by the (4).
where ∆y and ∆z are the cell dimensions; Z zy the impedance of the arms, equal to Z 0 = µ 0 /ǫ 0 for the standard SCN node; and Z s is the surface impedance (SI). In the case of a good conductor, an approximation by a real number for Γ is presented in [28]. However, Γ is in general complex and would alter the shape of the excitation pulses, which cannot be accounted for in the TLM method [7].

Application to confined environments modeled by arbitrary cross-section oversized waveguides: 2.5D TLM approach
The analysis of electrically large structures like tunnels involves a high computation time by using full-wave methods like TLM. However, some astutenesses can be considered to efficiently model them, which translates in finding an approach to simplify the analysis of oversized waveguides. From the geometrical point of view, tunnels can have any arbitrary shape, and this constraint has to be considered. From the electromagnetic point of view, waves propagating in tunnels can be classified in modes. This problem is well known in electromagnetic theory. As we shall see, a full-wave model based on a modal decomposition of the waves for the analysis of the radio propagation in tunnels at high frequencies is presented. Besides, the concept is used to further reduce the complexity of this electromagnetic problem.

The TLM modal approach
Modal approaches in frequency domain are based on the expansion of fields in terms of modes. In general, total fields can be represented as the sum of modes as shown in (5): where N is the total number of modes, A i (x, z) the amplitude of the i-th mode at the coordinate (x, z) in the plane perpendicular to the propagation direction y, α i its attenuation, and β i its phase constant. The mode amplitude A i (x, z) gives the electric or magnetic field distribution of the mode. The phase constant β represents the change in phase along the path traveled by the wave. Lastly, the attenuation constant is a measure of losses in the structure.
In the literature, the simplification of uniform guide problems in TLM was first introduced by [16] to obtain the dispersion properties. They proposed the reduction of the calculation region by introducing the field dependence solution in the propagation direction y. The longitudinal dependence can be described by exponential terms e −jβy , where β is the phase constant and y the longitudinal distance in the guide. Hence, for a specific mode, two points along the longitudinal distance of the guide have only a phase difference β (y 2 − y 1 ). Then, it is possible to find a relationship between the reflected voltages of the node at the time k∆t and the incident ones at the time (k + 1)∆t. This approach allows one to reduce the computation domain from a 3D mesh of 3D nodes to a 2D one, avoiding to calculate the fields over all points along the propagation direction.
The concept of mode will be employed to find a more efficient formulation, to obtain a reduced 2.5-dimensional TLM modal approach for finding fields. The term 2.5D is used as the 3D computational domain is reduced to a 2D mesh in the guide cross section. However, cells are 3D ones that account for all 6 electromagnetic field components. The reduced model can be applied for a uniform (invariant cross section along y) tunnel of arbitrary cross section.

Formulation of the 2.5D TLM node for guiding structures
To obtain the reduced TLM formulation for guiding structures, the voltages in each arm should be updated at each time step. In doing this, eight relationships between the reflected and incident voltages are needed. As mentioned before, stubs are added to model the background media with constitutive parameters ε and µ. Then, six additional relationships are needed for these stubs.
In [25], a procedure to find the required updating relationships based on Maxwell's equations is presented. To employ this procedure, six additional relationships are required to update all field components. The process is divided into two steps. The first step consists of finding the relationship between fields at the center of the node and the incident and reflected voltages taken at each arm. Then, the six relationships for the fields at instant n ∆t are found by using the integral form of Maxwell's equations in the three principal planes (z, y), (x, z), and (x, y). In the next step, field continuity equations at the center are applied by using the integral form of Maxwell's equations for contours passing through the center of the three principal planes. Then, 12 relationships for the voltages are found in terms of the fields at the center. By manipulating these expressions, the required eight relationships to update the voltages are found. Finally, the voltages for the stubs are calculated in the same way as the standard SCN node. The TLM algorithm is then completed.

Starting point: maxwell's equations in integral form
The starting point is the integral form of Maxwell's equations for an isotropic medium and given by (6) and (7).
where J e and J m are the electric and magnetic current density sources, respectively; the tensorsσ e ,μ,σ e ,σ m are given byε with ε 0 and µ 0 the permittivity and permeability of free space, respectively.

Field sampling
Consider the standard 3D SCN node. Using the same notation for the SCN node proposed by Johns, samples for the electric and magnetic fields in the cell are presented in (Figure 2). E n and H n correspond to field samples at each point of the cell. The incident V i and reflected V r voltages are related with these samples by (8). This equation gives us the relationship between the tangential fields at each face of the node and voltages in TLM network.

Determination of the updating equations for the fields
To find the updating relationships for the electric and magnetic fields, an approximation of Maxwell's equations in the integral form must be done. This can be accomplished by considering the rectangular mesh of Figure 2 and applying (6) and (7) in the main integration planes, (z, y), (x, z), and (x, y).  Maxwell-Ampère and Maxwell-Faraday equations applied to the contour Γ bounding this plane are given by (13) and (14), respectively. Moreover, electric and magnetic laws are given by (15) and (16), respectively. The calculation of these integrals is done by approximating their values as the mean value of the same integral at the time steps (n + 1/2)∆t and (n − 1/2)∆t. For the contour integrals, the approximation

Integration in the (z, y)-plane
Regarding the surface integrals, they can be calculated using where The application of the above expressions gives T corresponds to a delay operator, Z 0 is the characteristic impedance in free space, and c 0 the speed of light in vacuum. The termsẼ x andH x correspond to an approximation of the fields at instant n∆t. They can be calculated as the mean value of the field components in the plane:Ẽ whereŶ sx andẐ sx are the normalized admittances and impedances, respectively, to model the materials with a relative permeability and/or permittivity higher than unity. Their values are given byŶ Advanced Electromagnetic Waves 248 Regarding loss terms, the approximation of the corresponding integrals in (13) is done by taking their value at the center of the cell. This leads to where ζ denotes the field component x, y, or z and G eζ and R mζ correspond to the electric conductance and magnetic resistance, respectively, for each direction. Finally, the surface integrals of the sources are evaluated by taking their values at the center: where V eζ = J eζ ∆S and V mζ = J mζ ∆S. The infinitesimal area can be calculated, for instance, as ∆S = ∆y∆z for the x-direction (ζ = x).
By replacing the above expression in Maxwell's equations (6) and (7), the expressions in (24) and (25) can be obtained. Equation (8) was employed to express the fields in terms of reflected and incident voltages.

Integration in the (x, z)-plane
This plane is perpendicular to the propagation direction; thus, any integral is affected by the exponential dependence of the fields. They can be calculated as for the standard SCN node. The contours of integration for the magnetic and electric fields in the (x, z)-plane are shown in Figure 4, respectively.

Integration in the (x, y)-plane
This plane is parallel to the propagation direction. Thus, the value of the integrals will depend on the exponential dependence of the fields. The contours of integration for the magnetic and electric fields in the (x, y)-plane are shown in Figure 5.
In the above, G ez = σ ez ∆x∆y/∆z is the term used to model the electric losses, R mz = σ m ∆x∆y/∆z models the magnetic losses, V ez = J ez ∆x∆y, V my = J my ∆x∆y, and the terms E z andH z which correspond to the mean value of all field components in the plane are given, respectively, byẼ The normalized admittancesŶ sz and impedancesẐ sz are given by (19). By following a similar procedure to the other integration planes, the following expressions are obtained: The relationships of (24) z , from which fields can be computed at any time step. Eight additional relationships are still needed to obtain reflected voltages. By expressing the fields at the center of the cell, these expressions can be found, as detailed in the next subsection.

Determination of the updating equations for the voltages
To obtain the reflected voltages at any time step, six new contours have to be defined. Figure 6 shows the integration contours for finding fields at the center of the cell using Maxwell's equations in the integral form in the three planes.
The procedure is the same as presented before. First, Maxwell's equations are approximated for each of these contours. Consider the contour in dashed lines for the plane (z, y); Maxwell-Ampere and the electric flux equations are given, respectively, by Then, the expressions in (43) and (44) are used in combination with (8) to express the result in terms of reflected and incident voltages, yielding (45). By considering the remaining red and blue contours and applying Maxwell's equations in these planes, (46) to (50) can be found.
In summary, 12 expressions at the center were found. These expressions relate fields at instant n∆t, incident and reflected voltages. They allow us to update the eight voltages in the node. By replacing the above expressions of the fields in (8), one obtains The calculation of reflected voltages at stubs is done as for standard SCN:

Results of the 2.5D modal approach
To characterize any waveguide field amplitudes, attenuation constant and phase constant must be calculated. Here, we examine how to obtain these parameters by using the 2.5D TLM approach for guiding structures.

Phase constant: dispersion diagram β − ω
For each mode, β is a function of frequency ω; it is known as the dispersion diagram of the structure. In general, this relationship is not linearly proportional to the frequency for waveguides. The computation of the dispersion diagram with the TLM approach for guiding Advanced Electromagnetic Waves 254 structures is as follows: The structure has to be excited as explained in Section 3.1. Then, the field components at the time step k∆t can be found through (17,18,30,31,39,40). After a sufficient number of time-step iterations, the response corresponds to the superposition of the modal fields in the guide. By performing a Fourier transform, the frequency response of the structure can be obtained. The peaks in this spectrum correspond to the modes, and they belong to a specific value of the phase constant β. Hence, different β values have to be imposed to obtain the dispersion diagram.

Attenuation constant
The attenuation of the modes can be extracted from the time-domain response at a phase constant β l given by where M is the number of modes propagating in the structure, E lm (0) is the initial field magnitude of the m-th mode, Q lm its quality factor, and ω lm the angular resonance frequency of the mode in radians per second for the l-th evaluated value of β l . A similar expression can be found for the magnetic field by replacing E by H in (53). To characterize the attenuation of the modes, the damping parameter has to be estimated. At each frequency ω lm in the dispersion curve of a given mode, the attenuation constant is related with the quality factor Q lm and the group velocity V g lm by [34] The terms ω lm /2Q lm can be extracted from the sampled time-domain response of (53) by using an estimation technique for modal content of a time-varying waveform. For instance, the Matrix pencil method has been a very popular technique for efficient extraction of the modal parameters [14]. Finally, the group velocity in (54), which involves the derivate of β with respect to ω, can be approximated by taking the finite differences ∆β/∆ω of the dispersion diagram calculated in the previous section.

Implementation of the 2.5D TLM for guiding structures
To validate this approach, the study of a simple theoretical "reference solution" for a canonical geometry is considered. A metallic-rectangular waveguide homogeneously filled with a lossy dielectric material is studied in this section.

Results for a homogeneously lossy dielectric-filled metallic-rectangular waveguide
Consider a metallic waveguide with dimensions a = 6 cm and height b = 4 cm. It was filled with a lossy dielectric with relative permittivity ǫ r = 2 and a conductivity σ = 0.001 Sm −1 and simulated with the 2.5D TLM, as shown in Figure 7 (a). An array of 14 × 10 nodes is ✶✽ ❆❞✈❛♥❝❡❞ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡s considered, and a perfect electric boundary condition is imposed at the top, bottom, and side walls, as illustrated in Figure 7 (b). To reduce the dispersion error, a mesh size of ∆l = λ/15 at 6 GHz is chosen. The signal is obtained after N = 4000 time steps and is convolved with a Hanning window. The dispersion diagram is obtained in the frequency range from 1 GHz up to 4 GHz and is compared with the theory. For the case of a rectangular-metallic waveguide, the theoretical expression for the dispersion diagram is given by [26]: where ω is the angular frequency; µ = µ 0 ; ǫ = ǫ r ǫ 0 ; m, n correspond to the order of the modes; and a and b are the waveguide dimensions. The comparison between the numerical dispersion diagram and the theoretical one is shown in Figure 8.
The continuous lines represent the theoretical results, whereas the discrete points are the results obtained by using the 2.5D TLM algorithm. Different β values ranging from 0 to 80 (rad/m) with 80 samples have been considered. The simulation is carried out for 4000 iterations for each β value. Good agreement between theory and the 2.5D TLM for guiding structures was obtained. The dispersion curves of a plane wave propagating in free space and a wave propagating in a medium with the same relative permittivity as the loading material (ǫ r = 2) are plotted as well. It is worth observing that the dispersion curves of the modes tend to the curve of a plane wave propagating in a medium with relative permittivity ǫ r = 2 as the frequency increases. Finally, the theoretical attenuation constant of the metallic waveguide is calculated by [26]: where k = ω √ µǫ is the (real) wave number in the absence of losses, tanδ = (ωǫ ′′ + σ) /ωǫ ′ is the loss tangent with ǫ = ǫ ′ − jǫ ′′ = ǫ r ǫ 0 − jσ/ω, and β is the phase constant. The results Advanced Electromagnetic Waves 256  obtained for the attenuation constants by numerical simulation and theory are shown in Figure 9.
The comparison between both curves validates this approach quite satisfactorily. The discrepancies observed close to the cutoff frequency are explained by the fact that (57) is valid as long as α << β.
Consequently, since α is deduced from ∆β in the numerical approach, an inevitable large calculation error on the attenuation constant α is obtained for a given error on ∆β. To

Time-Domain Electromagnetic Wave Propagation in Confined Environments 257
✷✵ ❆❞✈❛♥❝❡❞ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡s illustrate this point, let us consider the calculation of the propagated error on α. The total differential of a function z = f (x 1 , x 2 , ..., x N ), where x i is the i-th independent variable, is given by the term dz denotes the error in z. The uncertainty in the calculation of the attenuation constant for waveguides can be found by considering the differential of the attenuation and phase constants, α and β, respectively. The differential of the attenuation and phase constants can be computed by Finally, by dividing (59) by (60), the differential of the attenuation constant in terms of the differential of the phase constant is given by is obtained.
As frequency increases, the first partial derivative on the right-hand side of (61) tends to zero and the second to the speed of light in free space, c 0 . Hence, for high frequencies, the error ∆α is bounded. For low frequencies, the first derivate tends to infinity and the second to a finite value, resulting in a larger error. Indeed, in Figure 9, the highest difference between both approaches occurs at cutoff frequencies of the modes and decreases far from these frequencies. Apart from this relative error on α, the comparison between numerical results and theory for both propagation constant β and the attenuation constant α confirms the suitability of the 2.5D TLM to characterize modal parameters in waveguiding structures. impedance boundary conditions (SIBCs), constitute the mathematical artifice to limit the computational domain.

Assumptions for tunnel wall modeling: SIBC concept
To apply the SIBC concept, fields in the walls region are assumed to be time exponential decaying, source free, and with no reflections returning to the interface air and tunnels walls. For practical purposes, the wall thickness t should be larger than the skin depth δ, (t >> δ). A thickness of 2δ would ensure an amplitude of 1.8% of the launched wave [30]. Another important requirement is the smallest radius of curvature on a surface R min , which must be larger than the skin depth (R min >> δ). Consequently, the interface can be considered locally flat and plane-wave concept can be applied. Considering a numerical method like TLM, using the common criteria of λ/10 for the maximum mesh size, waves can be treated locally as plane waves. Moreover, at this resolution level, the staircase mesh is a good approximation of curved surfaces.

SIBC implementation in time domain for tunnel wall modeling
The implementation of the SIBC concept was first introduced for good conductors. An overview of the implementations of these techniques in time domain can be found in [5,23,27]. Tunnel walls are typically modeled by a lossy dielectric with relative permittivity ǫ r = 10 − 12 and conductivity σ = 0.01 − 0.02Sm −1 [9]. Maloney in [21] first introduced the SIBC for lossy dielectric in FDTD by using an expansion of the impedance model proposed by [5] for good conductors. Finally, [27] proposed a way to improve its implement in FDTD.
In contrast with good conductors, the SIBC for lossy dielectric varies with the angle of incidence and polarization of incident waves. If one considers Figure 10, surface impedances for vertical and horizontal polarizations are given by [24] Time-Domain Electromagnetic Wave Propagation in Confined Environments 259 ✷✷ ❆❞✈❛♥❝❡❞ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡s Figure 10. Geometry and coordinate system for a TE and TM polarized waves incident onto a dispersive half space Z (s) is the surface impedance; s is the complex frequency; η is the intrinsic impedance of the lossy dielectric; ǫ, µ, and σ its constitutive parameters; and θ i and θ t are the angles of incidence and refraction, respectively. In tunnel environments, the incident angle θ i for the waves is unknown. Moreover, multiple reflections of waves with different angles of incidence are involved. To the best of the authors' knowledge, no work has been reported for SIBC under these conditions. A procedure to implement it in any time-domain method as well as its validation is presented in [2,4]. Here, a brief explanation of its implementation is described.

SIBC formulation for tunnel wall modeling
The SIBC is inherently a frequency-domain concept. When it is transformed into the time domain for use in methods such as TLM, it is replaced by a convolution integral. Maloney in [21] proposed a recursive formula to overcome the computational difficulties for its calculation. To implement the SIBC in TLM, the reflection coefficient of (4) has to be calculated. This coefficient depends on the polarization, frequency, and angle of incidence. Thus, their influences have to be studied.
For tunnel walls, the relative permittivity ǫ r >> 1. Then, for these values in (64), the term cos (θ t ) in (Eqs.62, 63) tends to one, obtaining the two angle-independent expressions of the surface impedance for both polarizations: By using a rational polynomial approximation [12], the resulting coefficient of (4) can be calculated for both polarizations into the form where ξ i and ζ i correspond to the poles and zeros ofP (s) and s i are the sample points.
Consider discrete values of a function P (S) at points S i . One can find two polynomials N (s) and D (s) of order M, such thatP (s) = N (s) /D (s) best approximates P (s), in the least square sense.
The next step is to obtain (67) in z-domain. The bilinear z-transform can be employed, obtaining the rational polynomial expression Finally, the discrete-time state variable representation of (68) can be obtained, which gives us the required relationship between the reflected and incident voltages necessary to implement the boundary in TLM. A canonical representation for the state-space and output equations is given, respectively, by This procedure guarantees a good agreement with theory for a relative permittivity ǫ r >> 1, i.e., it applies to roadway and railway tunnel modeling.

Wave Propagation in a Rectangular Tunnel with Lossy Dielectric Walls Using the SIBC
The wave propagation in a rectangular tunnel with transverse dimension w = 7 m × h = 5 m was studied. This tunnel was simulated using the TLM for guiding structures and time-domain SIBC implementations presented previously. The tunnel walls were modeled by using a lossy dielectric with ǫ r = 12 and σ = 0.02 Sm −1 . The tunnel is shown in Figure 11.
In a rectangular hollow dielectric waveguide, the modes are classified into two families: the first ones with most of its electric field polarized in the y-direction, designed as E y nm , and the other ones with most of its electric field polarized in the x-direction, designed as E x nm .  between our procedure and the commonly used Marcatilli's theory for dielectric rectangular waveguides (see Section 2.2). Some very good agreement was obtained for both modes. ) E x 11 mode of the rectangular tunnel calculated in the frequency range from 0.4 GHz to 2 GHz with Marcatilli's theory and our procedure, b) dispersion diagram for a rectangular tunnel in the frequency range from 0.2 GHz to 2 GHz using our numerical procedure [4] As it can be seen, the dispersion curves for the modes tend to the limit β = ω/c, and up to 3 GHz, they are propagating modes. They propagate with different phase constants and the propagation is strongly multimodal, e.g., around 120 modes at 3 GHz. The figure on the right of Figure 12 shows the comparison between both approaches in terms of attenuation constant. One can notice some discrepancies in the low region of the spectrum. One should stress that (54) is valid for α << β, and the approximation fails close to the cutoff frequency of the modes (see Section 4.1). Moreover, discrepancies can be explained by the fact that Marcatilli's theory neglects fields at the corners. As the frequency increases, fields are more concentrated in the core region, reducing this error. In Figure 12, the difference between the attenuation constants diminishes as the frequency increases, confirming this fact. Figure 13 shows the field profiles for the E y 11 and E x 11 modes in the waveguide calculated with Marcatilli's theory. Results obtained in Figure 14 were with the 2.5D TLM approach for guiding structures at 1 GHz. The structure was excited with a Gaussian pulse and injecting Advanced Electromagnetic Waves 262 ❚✐♠❡✲❉♦♠❛✐♥ ❊•❡❝tr♦♠❛❣♥❡t✐❝ ❲❛✈❡ Pr♦♣❛❣❛t✐♦♥ ✐♥ ❈♦♥✜♥❡❞ ❊♥✈✐r♦♥♠❡♥ts ✷✺ ✶✵✳✺✼✼✷✴✻✶✸✹✵ the corresponding phase constant of the mode, β ≈ 20 in Figure 12. Good agreement was found for the field profiles of both modes. and E x 11 modes calculated with Marcatilli's theory [18,19]. and E x 11 modes calculated with our procedure [4].
To conclude, it is important to mention that the application of this method can be extended to any uniform waveguiding structure, and also, it can be used to optimally determine antenna field specifications and positioning in confined or waveguide environments as presented in [3].

Conclusion
The TLM (transmission-line matrix) to characterize modes in guided structures was proposed and described. The 2.5D TLM approach is useful for finding the dispersion diagram, amplitude, and attenuation constant of modes propagating in uniform lossy guiding structures with arbitrary cross section. Its validation in a metallic-rectangular waveguide was presented. This confirms that this approach is well suited for numerical implementation.
Then, the modeling of the wave propagation in lossy waveguides as mine, roadway, or railway tunnels was studied. A procedure to implement the concept of SIBC (surface impedance boundary condition) to simulate the lossy dielectric walls was described. This concept allows us to avoid the meshing of the tunnel walls by replacing the wall medium by a frequency-dependent reflection coefficient. Thus, field values beyond the air-wall interface are not required for mode parameter computation. However, the presence of the dielectric beyond the interface is accounted for by the SIBC. The SIBC concept yields satisfactory results as long as ε r >> 1, for arbitrary conductivity. The illustration of the capabilities of this method in a rectangular tunnel shows that our procedure is in agreement with theory for the mode amplitudes, attenuation, and phase constants for the E y 11 and E x 11 modes when compared with Marcatilli's model. The 2.5D TLM approach for guiding structures and the SIBC concept might be also applied to study the radio-wave propagation for any cross section. The results in terms of amplitudes and dispersion curves and attenuation constants will allow to study different types of structures and the sensibility of these parameters with the polarization, excitation, shape, size, or frequency.