Finite-Region Approximation of EM Fields in Layered Biaxial Anisotropic Media

: A new algorithm is developed to accurately compute the electromagnetic (EM) ﬁelds in the layered biaxial anisotropic media. We enclose the computational region in an inﬁnitely long rectangular region by four vertical truncation planes and establish the corresponding algorithm to approximate the EM ﬁelds in the entire space. The EM ﬁelds in this region are expanded as a two-dimensional (2-D) Fourier series of the transverse variables. By using the spectral state variable method, the generalized reﬂection coefﬁcient matrices and transmission matrices are then derived to determine the Fourier coefﬁcients per layer. Therefore, we can obtain the spatial-domain EM ﬁelds by summing the 2-D Fourier series. To enhance the accuracy and efﬁciency of this algorithm, we apply the method of images to estimate the inﬂuence of the artiﬁcial boundaries on the EM ﬁelds at the observer. We then further develop a quantitative principle to choose the proper size of the region according to the desired error tolerance. With the proper choice, the summation of the series can achieve satisfactory accuracy. This algorithm is ﬁnally applied to simulate the responses of the triaxial logging tool in transversely isotropic and biaxial anisotropic media and is veriﬁed through comparisons to the other method.


Introduction
Electrical engineering and geophysical exploration applications involve electromagnetic (EM) fields in the layered (horizontally or cylindrically) media, including the design of microstrip circuits and antennas, airborne electromagnetic surveys, well logging, and so on [1][2][3][4][5][6][7]. In addition to the numerical methods such as the finite-difference and finitevolume method [8,9], the traditional integral-based analytic methods can also deal with the layered structures. For the horizontally layered media, one can express the spatial-domain EM fields as a two-dimensional (2-D) Fourier integral involving mixed spectral-domain EM fields (i.e., these fields are functions of one spatial variable and two spectral variables). When the medium is isotropic or transversely isotropic (TI), the spectral-domain EM field can be decomposed into TE and TM waves and can be further solved by introducing the scalar generalized reflection coefficients [10][11][12]. Then, one can obtain the spatial-domain EM fields by evaluating the so-called Sommerfeld integrals [13], which are derived from the 2-D Fourier integrals due to the cylindrical symmetry. However, when the media are biaxial anisotropic (BA), this decomposition will be invalid because the TE and TM waves are coupled to each other at interfaces.
The spectral state variable method has been developed to deal with the BA cases [14][15][16]. By defining a spectral state variable vector (which can be chosen differently as in [17,18]), one can derive a first-order differential system with a 4 × 4 system matrix. This system governs the spectral-domain EM fields, and its solutions for the homogeneous media can be solved easily. The solution for the layered media can be obtained by combining each homogeneous layer's solution with the generalized reflection coefficient matrices. Furthermore, one can also apply the 3-D Fourier transform to express the spatial EM fields as 3-D Fourier integrals involving the spectral-domain EM fields with a dependence of full spectral variables [19]. No matter which representation is used, the 2-D Fourier integrals always arise when one implements the spectral-to-spatial transformation. Unlike the integral for the TI cases, which for the BA cases cannot degrade into 1-D integral, its efficient evaluation is more challenging because of its high dimensionality. Sainath et al. [20] applied a complex-plane Gauss-Laguerre quadrature to evaluate this integral. To make the integrand decay fast, the singularity subtraction was introduced by Hu et al. [21], and a variation of this technique was provided by Hong et al. [18].
In recent decades, the problem of exploration and development of the anisotropic reservoirs, which are typical lossy media, has attracted a lot of attention in geophysical EM well logging [22][23][24][25][26] and the inversion of logging data [27,28]. To simplify the EM modeling in such media, the approximate solution of the EM fields is assumed to be a 2-D Fourier series of transverse variables x and y. This expression requires the entire space to be truncated by four vertical planes for the Fourier series and is only applicable to a finite interval. Therefore, the computational region is an infinitely long (along z-axis) rectangular region with a finite cross-section. In our previous work [29], we developed this finite-region approximation (FRA) technique to model the EM fields in the simple homogeneous media. In this paper, FRA are applied to deal with the layered BA media by combining it with the spectral state variable method. Furthermore, a quantitative principle for choosing the proper size of the region is introduced.
First, both the EM fields in this region and the source quantity are expanded as 2-D Fourier series about transverse variables x and y. The Fourier coefficients can be regarded as the discrete spectral-domain EM fields. Second, the mentioned spectral state variable method is then employed to determine the spectral-domain EM fields per layer. Finally, we can obtain the spatial-domain EM fields by summing the 2-D Fourier series. The error of the proposed FRA relative to the exact solution of the original problem in the entire space is the reflections from the four artificial truncation planes. Since the amplitude of EM waves in lossy media will decay rapidly as it propagates, the reflections will be negligible as long as the region is large enough. However, an oversized region may lead to the slow convergence of the series. To enhance the efficiency of FRA, we apply the method of images to quantitively estimate the error of the EM field at the observer for the homogeneous isotropic case and further develop a quantitative principle to determine the proper size of the region according to the desired error tolerance. Once the region's size is chosen, the proper truncation order of the summation of series can be determined by using a simple criterion in the process of summation. With the proper choice of the region's size, the summation of the series can achieve a satisfactory accuracy with a relatively small truncation order. Note that the spectral-to-spatial transformation in our new method is implemented by the simple summation of the 2-D Fourier series, without resorting to the additional numerical quadrature algorithms as in these integral-based methods.
In the numerical results section, the proposed FRA is applied to simulate the responses of the triaxial logging tool in the layered TI and BA media. The singularity subtraction similar to [13,21] is used to address the challenge of slow convergence for the highly deviated well. The agreements between results obtained by our method and those by the transmission line method (TLM) [30] validate our algorithm.

Spectral State Equation and Its Solution in Homogeneous Media
As shown in Figure 1, an infinitely long region Ω = [−L, L] × [−L, L] × (−∞, ∞) is occupied by a J-layer anisotropic medium with interfaces z j (j = 1, . . . , J + 1). Region Ω has a finite square cross-section with length 2L in the xy-plane. The permittivity ε, permeability µ of the medium are assumed to have the value for vacuum, and the medium can be characterized by the conductivity tensors σ 1 , σ 2 , . . . , σ J . The tensor σ for an arbitrary layer is assumed to be symmetric: (1) Remote Sens. 2022, 14, 3836 3 of 21 finite square cross-section with length 2 L in the xy-plane. The permittivity  , permeability  of the medium are assumed to have the value for vacuum, and the medium can be characterized by the conductivity tensors 1 2 , ,..., J σ σ σ . The tensor σ for an arbitrary layer is assumed to be symmetric: . xx xy xz xy yy yz xz yz zz Figure 1. The infinitely long rectangular region with a square cross-section and a layered anisotropic medium.
where ( , , ) x y z   r represents the position of the observer.
To find the solution of (2), we first expand ( ) E r and ( ) H r in terms of 2-D Fourier series about variables , x y : where n k n L   and m k m L   , , m n   (  is the set of integer number). In analogy to those integral-based methods, the unknown Fourier coefficients nm E  and nm H  can be regarded as the EM fields in the (discrete) spectral domain, and the integer pair nm is the spectral index. The delta function can also be expressed as a 2-D Fourier series [29], i.e.,  Given a magnetic dipole source M = M xx + M yŷ + M zẑ (x,ŷ,ẑ are the unit vectors in the three coordinate directions, and the time dependence is assumed to be e −iωt ) located at r ' = (0, 0, z ) ∈ Ω, the EM fields generated by this source satisfy Maxwell's equations: where r = (x, y, z) ∈ Ω represents the position of the observer.
To find the solution of (2), we first expand E(r) and H(r) in terms of 2-D Fourier series about variables x, y: where k n = nπ/L and k m = mπ/L, m, n ∈ Z (Z is the set of integer number). In analogy to those integral-based methods, the unknown Fourier coefficientsẼ nm andH nm can be regarded as the EM fields in the (discrete) spectral domain, and the integer pair nm is the spectral index. The delta function can also be expressed as a 2-D Fourier series [29], i.e., By substituting (3), (4) into (2), we can derive the equations aboutẼ nm andH nm for each spectral index, i.e., ∇ ×Ẽ nm e i(k n x+k m y) = iωµ 0 [H nm + Mδ(z−z ) and where S nm = [iωµM y , −iωµM x , ik m M z , −ik n M z ] T /4L 2 is the source vector; A nm is the 4 × 4 system matrix that can be partitioned into . These submatrices are The system matrix A nm has the factorization with eigenvalues in diagonal matrix Λ nm and eigenvectors in columns of matrix L nm . The 2 × 2 diagonal submatrices Λ u nm and Λ d nm consist of eigenvalues with positive and negative imaginary parts, respectively. The matrices L nm and Λ nm can be obtained by using linear algebra libraries or by analytically solving the eigenequation of A nm [15]; therefore, they can be regarded as the known quantities. Substituting (9) into (6) yields where represents the mode-wave vector and the new source term, respectively. In order to facilitate the derivation, we can divide Σ nm into sub-vectors For the homogeneous media, the solution of (10) can be expressed as The upper sub-vector of w nm represents the upward mode-wave since Λ u nm consists of eigenvalues with positive imaginary part. Similarly, the lower sub-vector represents the downward mode-wave.

Solution of Mode-Waves in the Layered Media
Based on the solutions of mode-waves in the homogeneous media and the superposition principle, we can further find the solutions in the layered media. The subscript nm signifying the spectral index will be omitted for the sake of brevity, and the layer index j and s are used to signify an arbitrary layer and the source layer, respectively. To facilitate the derivation, we introduce the following notations:

Formal Solution of Mode-Waves in the Source Layer
We first consider the source layer. To simplify the analysis, we assume that the source is located in one of the middle layers, i.e., z s < z < z s+1 , 1 < s < J. These horizontal interfaces will reflect the mode-waves produced directly by the source, which is the modewaves in the homogeneous medium. Thus, the incident mode-waves in this layer has the same form as (12), i.e., The reflected mode-waves can be expressed as where U s and D s represent the amplitudes of reflected upward wave at z = z s+1 and reflected downward wave at z = z s , respectively. Conversely, by introducing the generalized reflection coefficient matricesR j,j+1 andR j,j−1 [14], we can express U s and D s in terms of each other, i.e., Solving for U s and D s from (16) yields Combining (17) with (15), we can obtain the reflected part w ref s . Thus, for the source layer, we may calculate the incident wave and the reflected wave separately and then add these two parts to obtain the total mode-waves. For the convenient derivation of the mode-waves in the source-free layers, we should express w s in terms of the amplitudes of the total upward and downward mode-waves, as follows: Here, A − s is the amplitude for the total upward mode-waves at z = z s , A + s is the amplitudes for the total downward mode-waves at z = z s+1 . They have the following expressions:

Formal Solution of Mode-Waves in the Source-Free Layers
For the rest of the source-free layers, by using the generalized reflection coefficients [14,25], we can express the mode-waves as Here, expression (20) applies for layers above the source, and A − j is the amplitude for the total upward mode-waves at z = z j , whereas expression (21) applies for layers below the source, and A + j is the amplitudes for the total downward mode-waves at z = z j+1 . As in (18), (20) and (21), we have expressed the mode-waves in all layers in terms of the reflection coefficients and the amplitudes of the total upward and downward mode-waves. These unknown coefficients and amplitudes can be determined by matching boundary conditions at interfaces.
Suppose that the transmission coefficients T j,j+1 , T j−1,j and the local reflection coefficients R j,j+1 , R j−1,j are obtained in advance (see Appendix A). According to the field continuity conditions across interfaces z j below z , and we can derive the following equations: Solving for A + j andR j−1,j from (22) yields the following recurrence formulas Similarly, for these interfaces z j above z (j = 1, 2, . . . , s), we can derive the corresponding recurrence formulas Furthermore, we can allowR because there is no reflected downward mode-wave in the top layer and no reflected upward mode-wave in the bottom layer. These recurrence Formulas (23), (24), together with (19), (25) give the generalized reflection coefficients per interface and the amplitudes of mode-waves per layer. Then, the solution of mode-waves w nm for an arbitrary layer can be determined by using (18), (20) and (21).

Spatial-Domain EM Fields
Using the solution of w nm and the relation b nm = L nm w nm , we can obtain the four horizontal EM field components and further obtain the rest of the vertical components by using (7) and (8). In summary, the spectral-domain EM fields can be expressed as where F E nm and F H nm are 3 × 4 matrices, as expressed below Substituting (26) into (3) yields the spatial-domain EM fields, i.e., It can be easily verified by direct substitution that solution (28) satisfies the following boundary conditions: Hence, the FRA solution can be regarded as the solution to the Maxwell's equations in the region Ω subject to the boundary conditions in (29).

Tensor Green's Function and the Choice of Region's Size and Truncation Order
In the previous subsections, we derived the spatial-domain EM fields E(r), H(r) produced by a single magnetic dipole source Mδ(r − r ' ). Allowing M =x,ŷ,ẑ successively, we can compute the corresponding magnetic fields H x , H y , H z due to the three mutually orthogonal unit sources. Their combination produces the spatial-domain magnetic Green's function, i.e., where H pq denotes the p-component of the magnetic field due to the unit source in the q-axis direction. Using (28), we can also write the Green's function as a Fourier series, where the 4 × 3 matrix W nm = w nm,x , w nm,y , w nm,z represents the combination of the mode-waves corresponding to the three sources.
In practice, we can calculate the finite series in (31) by taking its partial sum, i.e., where the superscripts L and N are used to signify the region's size and the truncation order, respectively. Using the notation of partial sum, the FRA solution in (31) can be expressed as a limit when N → ∞ , i.e., When L → ∞ , the FRA solution further becomes the exact solution of the original problem in the entire space, Therefore, we will use the unified notation with double superscripts to express the can then be divided into two parts, represents the reflections from the four truncation planes, which is caused by the finite L, which we call the spatial truncation error, whereas is the spectral truncation error caused by the finite N. We first study the spatial truncation error.
For simplicity, we assume that the source is located at r ' 0 = (0, 0, 0) in an unbounded medium. As shown in Figure 2a is the solution to a free-space problem, whereas can be regarded as the solution to a boundary-value problem in Ω subject to the boundary conditions (29). Using the method of images, we can construct an equivalent free-space problem to the problem in Figure 2b by introducing an infinite number of image sources with the same magnitude and direction as the original source. As shown in Figure 2c, these sources are located at The one located at 0,0 0      r r represents the original source (shown in red dot) and the rest located at pq    r represent the image sources (shown in green dots). It is easy to verify that the sum of the contributions due to the original source and its image sources satisfies the conditions in (29). By the uniqueness theorem, we conclude that the problem in Figure 2c is equivalent to the problem in Figure 2b.
To measure the magnitude of a tensor, we introduce the following Frobenius norm for an arbitrary 3 3  tensor A with entries , , 1, 2,3 ij a i j  : Taking the norm of (38) and using the triangle inequality, we have The right-hand side of (40) gives an upper bound for the spatial truncation error. When the medium is isotropic, each term of can be calculated by using the exact formula of the magnetic Green's function in the free space [31]: r r , ˆˆˆˆˆ   I xx yy zz is the identity dyad. We use (41) to estimate the influence the computational region's size L on the EM fields and further develop a principle to choose the proper L even for the anisotropic media. The one located at r 0,0 = r 0 ∈ Ω represents the original source (shown in red dot) and the rest located at r pq / ∈ Ω represent the image sources (shown in green dots). It is easy to verify that the sum of the contributions due to the original source and its image sources satisfies the conditions in (29). By the uniqueness theorem, we conclude that the problem in Figure 2c is equivalent to the problem in Figure 2b. Thus, the FRA solution G L,∞ can be expressed as From (37), we obtain the following expression for the spatial truncation error: To measure the magnitude of a tensor, we introduce the following Frobenius norm for an arbitrary 3 × 3 tensor A with entries a ij , i, j = 1, 2, 3: Taking the norm of (38) and using the triangle inequality, we have The right-hand side of (40) gives an upper bound for the spatial truncation error. When the medium is isotropic, each term of can be calculated by using the exact formula of the magnetic Green's function in the free space [31]: where k = iωµ(σ − iωε), R = r − r ' , I =xx +ŷŷ +ẑẑ is the identity dyad. We use (41) to estimate the influence the computational region's size L on the EM fields and further develop a principle to choose the proper L even for the anisotropic media. Consider a specific case as an example: the observer is located at r = 1.016 × (1/2, 0, √ 3/2) m, and the frequency of the source is 20 kHz. The right-side of (40) only depends on the variables σ and L; thus, we can define it as the following function: Its P-th (P = 1, 2, . . . , ∞) order approximation can be defined as which represents the contributions of the (2P + 1) 2 − 1 image sources nearest to the original source. Figure 3 shows the graph of f P (σ, L) when P is equal to four different values 1, 2, 3 and 1000 in the case of σ = 1 S/m. From the results, we can observe that f P (σ, L) converges quickly with the increase in P and the values of f 3 (σ, L) becoming almost consistent with that of f 1000 (σ, L). Thus, we choose f 3 (σ, L) to approximate f (σ, L) to investigate the influence of region size L. The whole graph of f 3 (σ, L) in Figure 4 demonstrates that f 3 (σ, L) decreases monotonically with the decrease in σ and L. If the conductivity σ is known and some error tolerance denoted by e tol is given, we can find a unique proper value of L, denoted by L p , such that       We can observe that there exists a proper truncation order N p for each case of L p , and the total error of G L p ,N no longer decreases after N is up to N p . This observation indicates the following results: That is, G L p ,N p has arrived at the value of G L p ,∞ , and the spectral truncation error can be ignored. Consequently, the total error of G L p ,N p is equal to the spatial truncation error, which will not exceed the present error tolerance e tol , i.e., The numerical results in Figure 5 confirm this conjecture. In a practical calculation, we can determine N p as the smallest N that meets the following stopping criterion: and take G to the anisotropic case as well. Considering a TI medium of σ = diag(10 , 10, 1) S/m, Figure 6 presents the total error of G L p ,N as the truncation order N increases for the three L p determined by e tol = 10 −4 , 10 −5 , 10 −6 . The total error is calculated by where G ∞,∞ TI is the exact solution of the magnetic Green's function in an unbounded TI medium, and its formula can be found in [32]. For this medium, although we cannot guarantee that the error of approximate solution G L p ,N p is less than the present e tol , the results in Figure 6 demonstrate that the error is close to the present e tol . Thus, we still apply previous principles to determine L p and N p for the anisotropic media. For a J-layer anisotropic medium, we may take min σ 1 , . . . , σ j , . . . , σ J as the approximate isotropic conductivity, where σ j is the minimum of the three principal components of the conductivity in layer j. In the following section, the error tolerance will be fixed at e tol = 10 −5 .
is the exact solution of the magnetic Green's function in an unbounded TI medium, and its formula can be found in [32]. For this medium, although we cannot guarantee that the error of approximate solution p p , L N G is less than the present tol e , the results in Figure 6 demonstrate that the error is close to the present tol e . Thus, we still apply previous principles to determine p L and p N for the anisotropic media. For a J-layer anisotropic medium, we may take 1 min{ ,..., ,..., } j J       as the approximate isotropic conductivity, where j  is the minimum of the three principal components of the conductivity in layer j . In the following section, the error tolerance will be fixed at

Results
In this section, we apply the proposed FRA to simulate the responses of the triaxial logging tool in the layered TI and BA formations. As shown in Figure 7a, the tool is equipped with three mutually orthogonal transmitters T , T , T and three mutually or-

Results
In this section, we apply the proposed FRA to simulate the responses of the triaxial logging tool in the layered TI and BA formations. As shown in Figure 7a, the tool is equipped with three mutually orthogonal transmitters T x , T y , T z and three mutually orthogonal receivers R x , R y , R z . Figure 7b illustrates the orientation of the tool coordinates x t y t z t with respect to the formation coordinates xyz: the y t -axis is assumed to be parallel to the xy-plane and the z t -axis points to the direction from the transmitters to receivers. The rotation matrix from xyz to x t y t z t can be expressed as where α and θ represent the azimuthal and dip angles of the tool. To simulate the responses of the triaxial logging tool, the magnetic Green's function in (30) should be converted into the tool coordinates by H will be omitted for brevity. As shown in Figure 8, a five-layer anisotropic model is assumed, where the thi nesses of the three middle layers are 2, 2 and 4 m, respectively. The layers 1, 3 and 5 assumed to be isotropic, and the conductivities are 0.1, 0.1 and 0.05 S/m, respectively. T conductivities in layers 2 and 4 share the identical tensor conductivity σ , which is sumed to be TI or BA. As shown in Figure 8, a five-layer anisotropic model is assumed, where the thicknesses of the three middle layers are 2, 2 and 4 m, respectively. The layers 1, 3 and 5 are assumed to be isotropic, and the conductivities are 0.1, 0.1 and 0.05 S/m, respectively. The conductivities in layers 2 and 4 share the identical tensor conductivity σ, which is assumed to be TI or BA.

Triaxial Logging Responses in Layered TI Media
To validate this new algorithm, we compare the results by the present FRA with those of the TLM [30] in a TI formation. The five-layer model in Figure 8 is considered, where the conductivity σ of layers 2 and 4 is assumed to be TI and is equal to diag(1, 1, 0.1) S/m. Two different tool dip angles θ = 60 • and θ = 89 • are considered. We will only provide the results of the components ImH xx , ImH yy , ImH zz , ImH xz , ImH zx because the remaining ImH xy , ImH yx , ImH yz , ImH zy are all zeros. Figures 9 and 10 show the comparison with the TLM for the cases of θ = 60 • and θ = 89 • , respectively. The corresponding relative errors are displayed in Figures 9f and 10f. The relative error of each component is less than 2% except for a few sampling points (where the absolute value of response is close to zero). The excellent agreements validate our algorithm. The direct summation of the Fourier series for the case of θ = 89 • exhibits a weakly convergent behavior, and the required summation order will be extremely large. This difficulty is solved by using the singularity subtraction [13,21]. Furthermore, from Figures 9 and 10, we can find that the local extremums of the cross components ImH xz and ImH zx can clearly indicate the location of the interfaces for both cases.  Figures 9f and 10f. The relative error of each component is less than 2% except for a few sampling points (where the absolute value of response is close to zero). The excellent agreements validate our algorithm. The direct summation of the Fourier series for the case of 89    exhibits a weakly convergent behavior, and the required summation order will be extremely large. This difficulty is solved by using the singularity subtraction [13,21]. Furthermore, from Figures 9 and 10, we can find that the local extremums of the cross components Im xz H and Im zx H can clearly indicate the location of the interfaces for both cases.
To further validate the present algorithm in formation with a high conductivity, we assume the conductivity of layers 2 and 4 to be diag(20, 20,2) S/m  σ and the tool dip angle to be 60    . Figure 11 shows the comparison with the TLM. The agreement validates our algorithm again.    are considered. All of them share the same principal conductivity diag(4, 1, 0.5) S/m. The orientations of their principal coordinates x p y p z p with respect to the formation coordinates xyz are also illustrated in Figure 12: the principal coordinates x p y p z p of the simple BA case are consistent with xyz; those of the dipping BA case have an anisotropic dip angle γ = 15 • ; those of the Full-tensor BA case have an anisotropic azimuthal angle β = 15 • and an anisotropic dip angle γ = 15 • . The tool dip angle is fixed at θ = 60 • . For the simple BA and dipping BA cases, Figure 13 only displays the results of components ImH xx , ImH yy , ImH zz , ImH xz , ImH zx , because the rest are all zeros. Note that in these two cases, the tool dip angle and the anisotropic dip angle are in the same xz-plane. However, in the cases of full-tensor BA conductivity, this condition no longer holds because of the anisotropic azimuthal angle. The results in Figure 14 show that none of the nine responses are zero. Furthermore, comparing the results of the dipping BA (γ = 15 • ) case with those of the simple BA (γ = 0 • ) case, we can observe that the three main components ImH xx , ImH yy , and ImH zz are sensitive to the change in the anisotropic dip angle γ. These differences among the results of the three cases indicate that the orientation of BA conductivity has a significant influence on the responses of the triaxial logging tool. To obtain a more reliable interpretation of logging data, this effect should be considered.    Informed Consent Statement: Not applicable.

Data Availability Statement:
The data included in this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
By the proper normalization, the eigenvector matrix for layer j and its inverse matrix can be expressed in terms of the four 2 × 2 matrices L j,I , L j,II , L j,III , L j,VI (the detailed expressions can be found in [15]) and they transpose the matrices as follows: Consider a two-layer model and assume an incident downward wave with amplitude A + hitting the interface z 2 from upper layer 1 into lower layer 2. The mode-waves in these two layers are where R 1,2 and T 1,2 are the unknown local reflection and transmission coefficients, respectively. To determine these coefficients, we apply continuity conditions across the interface z 2 , i.e., b(z − 2 ) = b(z + 2 ) (A4) where z − 2 (z + 2 ) represent the left (right) limit of z 2 . Using the relation b = Lw and expressions in (A2) and (A3), the condition can be written as Similarly, assuming that an incident upward wave hits interface z 2 from the lower layer 2, and after a similar derivation, we obtain