**1. Introduction**

When wind waves generated in the deep-water approach coastal regions, they experience various physical phenomena caused by inferences with structures, bathymetric variations, nonlinear wave interactions, etc. To prevent coasts from huge wave attacks, either floating or submerged breakwaters are usually installed in the coastal area. Although the nonlinear effects become significant as the waves approach the shoreline, consistent linear solutions are still valuable and provide extensive information concerning the wave impact on the nearshore and coastal environments. Furthermore, the linear solution usually serves as the starting point for a weakly nonlinear model [1].

Both submerged and floating breakwaters are typically designed in coastal regions. Submerged breakwaters are conventional structures that often rest on the sea floor. In addition, they are built as submerged types to satisfy the requirements of coastal landscapes and ecologies. Submerged breakwaters are heavy and large, and are designed by engineers for different purposes [2]. On the other hand, floating breakwaters have the advantage of lower construction costs compared to submerged breakwaters. Floating structures may be used in fishing farms for ecology conservation, tourism, and leisure. Moreover, they may be implemented at the ocean engineering working stations, such as oil exploration stations [3]. Engineers also use assembled floating structures to construct an airport on the sea [4], or to provide hospitable environments on the surface of the water.

The Bragg reflection caused by the periodic breakwaters can help to effectively attenuate waves. For example, Mei et al. [5] considered using a series of submerged sinusoidal sandbars to protect the

ocean platforms in Ekofisk of the North Sea, and Bailard et al. [6] found that the Bragg reflection of submerged bars can protect U.S. Gulf Coasts from storm-induced waves. In addition, Tsai and Wen [7] indicated that the Bragg reflections of submerged breakwaters were e ffective for protecting the Mi-Tuo Coast, Taiwan. The vortex generation and dissipation accompanying the Bragg scattering of water waves propagating over a series of submerged rectangular breakwaters were investigated by Hsu*,* et al. [8]. Recently, the Bragg reflections of floating breakwaters were studied by Ouyang et al. [9] and Ding et al. [10]. In this study, the combined Bragg reflections of periodic submerged and floating breakwaters are considered.

Numerical solutions are inevitable for solving water wave-scattering problems as analytic solutions are rare [11]. Berkho ff [12] derived the mild-slope equation (MSE) by integrating the governing equation over the vertical interval of water depth. Subsequently, the MSE was modified and improved in various studies [13,14]. In addition to the prescribed one-equation models, further improvements were made by including the evanescent modes, resulting in a system of hierarchical MSEs [15,16]. Athanassoulis and Belibassakis [17] additionally included a sloping-bottom mode to formulate the consistent coupled-mode system (CCMS), which has been applied to many water wave problems [18–20]. The MSE has been applied to solve problems of nonlinear waves [18,21], three dimensions [22], wave–current interactions [21,23], time evolutions [24], Bragg reflections [13], floating structures [25,26], etc. A comprehensive review can be found in a recent article [27].

Alternately, Takano [28] developed the eigenfunction marching method (EMM) for solving normal incident wave scattering over an elevated sill and a fixed surface obstacle. Subsequently, Kirby and his coauthors [29,30] applied the EMM to solve problems of wave scattering over a trench of oblique incidences. For waves propagating over an arbitrary bottom profile, Devillard et al. [31], O'Hare and Davies [32,33], and Tsai et al. [34,35] decomposed the bottom profiles into a sequence of flat shelves separated by steps. The EMM has been applied to problems of viscous wave scattering [36–38], water wave scattering by tension-leg structures [39] and thin floating plates [40]. The accuracy of the EMM solutions was shown to be comparable with that of the MSE solutions [41]. In addition, the mathematical formulation is EMM is simpler as there are no requirements for the spatial derivatives of the eigenfunctions; however, they are needed in the MSE. However, the applications of EMM to three-dimensional, nonlinear, and/or time-dependent problems require further investigation.

In 1966, Kato, et al. [ ¯ 42] conducted laboratory experiments on the reflections of wave scattering using four simple floating structures, including a rectangular structure. Through numerical method, the di ffraction of oblique waves scattering by a surface-piercing rectangular structure was studied by Bai [43]. Sequentially, Kanoria et al. [44] derived analytical solutions for normally incident wave scattering by a surface-piercing rectangular structure in water of uniform finite depth. The analytical solutions were then extended to oblique waves by Söylemez and Gören [45]. For an arbitrary cross-section, Garrison [46] developed a Green's function procedure to compute oblique wave scattering. Using numerical methods, Ouyang et al. [9] recently studied the Bragg reflections of normal waves by fixed rectangular structures. Ding et al. [10] studied the Bragg reflections of normal waves by structures of di fferent shapes using the boundary element method. All the prescribed studies consider scattering problems with di fferent configurations of various structures over a flat bottom.

Manisha et al. [47] recently developed a model considering the e ffects of bottom undulations for oblique wave interaction with a surface-piercing rectangular structure behind a submerged breakwater or a trench. They connected the solutions of the MSE and EMM for the regions of the undulated bottom and rectangular structure over the flat bottom, respectively. In this study, the EMM model is developed for analyzing the combined phenomena of oblique incidence, surface-piercing structures of di fferent shapes, Bragg reflections, and undulated bottoms. In the solution procedure, the surface-piercing structures and bottom topography are sliced into successive flat shelves separated by abrupt steps. The matching conditions of the normal flow flux and the continuity of pressure are imposed on the interface boundaries. The EMM model is validated by comparisons with analytical solutions in the literature [10,45,47].

Bragg's law is usually used to predict the wavelengths at which the X-rays are intensively reflected by crystalline solids [48]. For water-wave problems, Bragg's law is applied for scatterings by floating [9,10] and submerged [5–7,23] structures at normal incidence. For oblique incidence, Mei [49] and Dalrymple et al. [50] applied Bragg's law for water wave scattering by bottom ripples. In this study, numerical experiments were conducted to study the Bragg reflections by the combined effects of floating structures, bottom variations, and oblique incidence. In addition, the numerical results are compared with those predicted by Bragg's law.

This paper is organized as follows: the wave problem is mathematically modeled and the EMM solution is developed in Section 2, and the EMM model is validated in Section 3. Discussions on oblique Bragg reflections by surface-piercing and submerged breakwaters are provided in Section 4. Finally, conclusions are presented in Section 5.

#### **2. Materials and Methods**

#### *2.1. The Mathematical Model*

We consider the problem of oblique monochromatic water wave scattering by surface-piercing structures over uneven bottoms. The wave amplitude is assumed to be small enough that the linear wave theory is applicable. The wave motion is assumed to be time-harmonic, *<sup>e</sup>*<sup>−</sup>**i**σ*t*, where σ = 2π/*T* is the angular frequency, *T* is the wave period, *t* is the time, and **i** is the unit of complex numbers. Figure 1 shows a schematic representation of the wave scattering problem induced by a surface-piercing structure over uneven bottoms. In the figure, the surface-piercing structures and sea bottom are discretized into a series of *M* shelves in the intervals of *xm*−1 ≤ *x* ≤ *xm* for *m* = 1, 2, 3 ... , *M*, with a water depth *hm* > *dm*, where *dm* > 0 is the submergence depth of the structure. Alternatively, *dm* = 0 is designated if there is no structure in the interval. Furthermore, *x*0= −∞ and *xM*= ∞ are assumed.

Considering the solution on the *m*-th shelf in the interval *xm*−1 ≤ *x* ≤ *xm*, the velocity of the fluid is defined by

$$\mathbf{u}\_m = \nabla \phi\_{m\prime} \tag{1}$$

where ∇ = (∂/∂*<sup>x</sup>*, ∂/∂*y*, ∂/∂*z*) is the three-dimensional del operator with respect to the three-dimensional Cartesian coordinates (*x*, *y*, *z*) and φ*m* is the velocity potential. According to the linear wave theory, the velocity potential is governed by the Laplace equation as

$$
\nabla^2 \phi\_m = 0.\tag{2}
$$

It should be noted that φ*m* is only related to the spatial part of the velocity potential for the remainder of this work. If there is no structure in the interval (*dm* = 0), the problem is subjected to the kinematic and dynamic free-surface boundary conditions, respectively, as

$$-\mathbf{i}\sigma\eta\_m - \frac{\partial\phi\_m}{\partial z} = 0\tag{3}$$

and

$$-\mathbf{i}\sigma\phi\_m + \mathbf{g}\eta\_m = 0 \text{ on } z = 0. \tag{4}$$

Equations (3) and (4) can be combined to obtain

$$\frac{\partial \phi\_m}{\partial z} - \frac{\sigma^2}{g} \phi\_m = 0 \text{ on } z = 0. \tag{5}$$

On the other hand, if there is a structure in the interval (*dm* > 0), the boundary condition on the bottom of the surface-piercing structure is given by

$$\frac{\partial \phi\_m}{\partial z} = 0 \text{ on } z = -d\_m. \tag{6}$$

In addition, the boundary condition on the sea bottom can be expressed in the form of

$$\frac{\partial \phi\_{\rm m}}{\partial z} = 0 \text{ on } z = -h\_{\rm m}. \tag{7}$$

**Figure 1.** Schematic representation of the boundary-value problem of water-wave-scattering by the surface-piercing structure over even bottom. The study domain divided into different regions with *M* shelves separated by *M* − 1 steps.

Equations (2)–(7) are sufficient to construct the complete solution by eigenfunctions and will be given in the next section.

Then, connection conditions are required to match the solutions φ*m* and φ*m*+<sup>1</sup> as

$$
\phi\_m = \phi\_{m+1} \tag{8}
$$

and

$$\frac{\partial \phi\_m}{\partial \mathbf{x}} = \frac{\partial \phi\_{m+1}}{\partial \mathbf{x}}\tag{9}$$

on −*d*max *m* < *z* < −*h*min *m* and *x* = *xm* with *d*max *m* = max(*dm*, *dm*+<sup>1</sup>) and *h*min *m* = min(*hm*, *hm*+<sup>1</sup>). In addition, no-penetration conditions are needed on the bottom and structure side-walls, respectively, as

$$\frac{\partial \phi}{\partial \mathbf{x}} = 0 \text{ on } -h\_m^{\text{max}} < z < -h\_m^{\text{min}} \text{ and } \mathbf{x} = \mathbf{x}\_{\mathcal{W}} \tag{10}$$

and

$$\frac{\partial \phi}{\partial \mathbf{x}} = 0 \text{ on } -d\_{\mathbf{m}}^{\text{max}} < z < -d\_{\mathbf{m}}^{\text{min}} \text{ and } \mathbf{x} = \mathbf{x}\_{\mathbf{m}}.\tag{11}$$

The definitions of *h*max *m* and *d*min *m* are similar, and thus neglected here. In Equations (10) and (11), φ stands for either φ*m* or φ*m*+<sup>1</sup> depending on the water side of the wall.

Then, considering a monochromatic wave train with incidence angle α, amplitude *a*, frequency σ, and wavelength λ, which propagates towards the surface-piercing structures over an uneven bottom. Therefore to make the solution unique, the following far-field conditions are required

$$\eta = a \Big(e^{\mathrm{i}\vec{k}\_{1,0}\cdot\mathbf{x}} + K\_R e^{\mathrm{i}\partial\_R} e^{-\mathrm{i}\vec{k}\_{1,0}\cdot\mathbf{x}}\Big) e^{\mathrm{i}\mathbf{k}\_y\cdot\mathbf{y}} \text{ as } \mathbf{x} \to -\infty \tag{12}$$

and

$$
\eta = aK\_T e^{\mathbf{i}\vartheta\_T} e^{\mathbf{i}\hat{\mathbb{k}}\_{M,0} \mathbf{x}} e^{\mathbf{i}\mathbf{k}\_y y} \text{as } \mathfrak{x} \to \infty. \tag{13}
$$

where *KR*, θ*R*, *KT*, and θ*T* are real numbers, such that *KRe***<sup>i</sup>**θ*R* and *KTe***<sup>i</sup>**θ*T* are the reflection and transmission coefficients, respectively. In Equations (12) and (13), ˆ *k*1,0, ˆ *kM*,0, and *ky* are positive real wavenumbers defined by

$$k\_{m,n} = \sqrt{k\_{m,n}^2 - k\_y^2} \,\text{.}\tag{14}$$

and

$$k\_{\mathcal{Y}} = k\_{\mathcal{I}, \mathcal{O}} \sin \alpha. \tag{15}$$

where *k*1,0 = 2π/λ > 0 and *kM*,<sup>0</sup> > 0 are the progressive wavenumbers obtained from the dispersion relation

$$\frac{\sigma^2}{\mathcal{S}} = k\_{m,0} \tanh k\_{m,0} l\_m. \tag{16}$$

Here, it is assumed that there is no structure over the first and last shelves, i.e., *d*1 = 0 and *dM* = 0.

Now, the problem is well-defined, and the estimation of the reflection and transmission coefficients are presented in the next subsection.

#### *2.2. Eigenfunction Matching Method*

To construct the complete solution using eigenfunctions, a complete set of wavenumbers is required. When there is no structure on a shelf (*dm* = 0), in addition to the progressive wavenumber *km*,<sup>0</sup> in Equation (16) the evanescent wavenumbers *km*,*<sup>n</sup>* (*n* = 1, 2, 3, ...) are defined by

$$k\_{m,n} = \mathbf{i} \kappa\_{m,n} \tag{17}$$

where <sup>κ</sup>*m*,*<sup>n</sup>* is the *n*-th smallest positive root of the dispersion relation

$$\frac{\sigma^2}{\mathcal{S}} = -\kappa\_{\mathfrak{m},\mathfrak{u}} \tan \, \kappa\_{\mathfrak{m},\mathfrak{u}} l\_{\mathfrak{m}}.\tag{18}$$

and when there is a structure over a shelf (*dm* > 0), the wavenumber is alternatively defined as

$$k\_{m,n} = \frac{\mathbf{i}n\pi}{h\_m - d\_m} \tag{19}$$

for *n* = 0, 1, 2, .... As the incident wave is oblique, we have to define the x-component of the wavenumber ˆ *km*,*<sup>n</sup>* by Equation (14). The complex-valued wavenumbers, *km*,*<sup>n</sup>* and ˆ *km*,*n*, defined in Equations (14)–(19), enable the formulation of a unified EMM for all types of situations (*dm* = 0 or *dm* > 0; *n* = 0 or n > 0).

Based on the linear wave theory, the complete solution of the velocity potential for the *m*-th shelf may be expressed as

$$\phi\_{\rm{fl}}(\mathbf{x}, y, z) = \sum\_{n=0}^{N} \left( A\_{m, \rm{n}} \boldsymbol{\zeta}\_{m, \rm{n}}^{(1)}(\mathbf{x}) + B\_{m, \rm{n}} \boldsymbol{\zeta}\_{m, \rm{n}}^{(2)}(\mathbf{x}) \right) \boldsymbol{\zeta}\_{m, \rm{n}}(\mathbf{z}) e^{i \mathbf{k}\_{g} y} \tag{20}$$

for *m* = 1, 2, 3, ... , *M*, where *Am*,*<sup>n</sup>* and *Bm*,*<sup>n</sup>* are unknown coefficients to be determined. Upon applying the conditions in Equations (2), (5)–(7) and by employing the method of the separation of variables, the eigenfunctions, ζ*<sup>m</sup>*,*<sup>n</sup>*(*z*), ξ(1) *<sup>m</sup>*,*<sup>n</sup>*(*x*), and ξ(2) *<sup>m</sup>*,*<sup>n</sup>*(*x*), can be obtained and expressed as

$$
\zeta\_{m,n}(z) = \cosh k\_{m,n}(h\_m + z),
\tag{21}
$$

$$
\xi\_{m,n}^{(1)}(x) = \begin{cases}
 e^{i\vec{k}\_{m,n}(x-\overline{x}\_{m-1})} & \hat{k}\_{m,n} \neq 0 \\
 1 & \hat{k}\_{m,n} = 0
\end{cases} \tag{22}
$$

and

$$\xi\_{m,n}^{(2)}(x) = \begin{cases} e^{-i\hat{k}\_{m,n}(x-\overline{x}\_m)} & \hat{k}\_{m,n} \neq 0\\ \infty & \hat{k}\_{m,n} = 0 \end{cases} \tag{23}$$

with

$$\begin{cases} \overline{\mathbf{x}}\_{\mathfrak{m}} = \mathbf{x}\_{\mathfrak{m}} \text{ for } \mathfrak{m} = 1, 2, \dots, M - 1\\ \overline{\mathbf{x}}\_{0} = \overline{\mathbf{x}}\_{M} = 0. \end{cases} \tag{24}$$

By observing at Equation (21), we have ζ*<sup>m</sup>*,*<sup>n</sup>* = 1 for *km*,*<sup>n</sup>* = 0. According to the Sturm–Liouville theory [51], the following orthogonal relation is used for solving the problem

$$
\left\langle \zeta\_{m,n} \middle| \zeta\_{m,l} \right\rangle = \Lambda\_{m,n} \delta\_{nl} \tag{25}
$$

where *n* and *l* is a mode index varying from 0 to *N*, δ*nl* is the Kronecker delta function, and Λ*n* is a function of *hm* and *km*,*n*, written as

$$
\Lambda\_{m,n} = \frac{2k\_{m,n}(h\_m - d\_m) + \text{sirk2}k\_{m,n}(h\_m - d\_m)}{4k\_{m,n}} \tag{26}
$$

For convenience, we define the inner product of two depth eigenfunctions as follows.

$$
\langle F|G \rangle = \int\_{-\lambda\_2}^{-\lambda\_1} F(z)G(z)dz,\tag{27}
$$

where *F* and *G* are the depth eigenfunctions of ζ*<sup>m</sup>*,*<sup>n</sup>* with arbitrary *m* and *n*; and λ1 and λ2 represent the structure submergence and water depths, respectively, which correspond to the first depth eigenfunction *F*.

It should be noted that the eigenfunction definitions of ζ*<sup>m</sup>*,*<sup>n</sup>*(*z*), ξ(1) *<sup>m</sup>*,*<sup>n</sup>*(*x*), and ξ(2) *<sup>m</sup>*,*<sup>n</sup>*(*x*) are valid for all cases (*dm* = 0 or *dm* > 0; *n* = 0 or *n* > 0) if the complex-valued wavenumbers, *km*,*<sup>n</sup>* and ˆ *km*,*n*, are defined by Equations (14)–(19).

Based on the far-field conditions (Equations (12) and (13)) and the dynamic free-surface boundary condition (Equation (4)), the far-field solutions of the velocity potential can be expressed as

$$\phi\_1 = -\frac{\mathrm{iag}}{\sigma} \frac{\cosh k\_{1,0}(h\_1 + z)}{\cosh k\_{1,0}h\_1} \Big( e^{\mathrm{i}\hat{k}\_{1,0}x} + K\_R e^{\mathrm{i}\partial\_R} e^{-\mathrm{i}\hat{k}\_{1,0}x} \Big) e^{\mathrm{i}k\_gy} \text{ as } x \to -\infty \tag{28}$$

and

$$\phi\_{\rm M} = -\frac{\mathrm{i}ag}{\sigma} \frac{\cosh k\_{\mathrm{M},0}(h\_{\mathrm{M}} + z)}{\cosh k\_{\mathrm{M},0}h\_{\mathrm{M}}} \Big(K\_T \epsilon^{\mathrm{i}\rho\_T} \epsilon^{\mathrm{i}\hat{k}\_{\mathrm{M},0}x} \Big) e^{\mathrm{i}k\_{\mathrm{j}}y} \text{ as } x \to \infty. \tag{29}$$

Comparing Equations (20), (28), and (29), the following equations can be obtained as

$$A\_{1,0} = -\frac{\text{igg}}{\sigma} \frac{1}{\cosh k\_{1,0} h\_1} \,\text{}\tag{30}$$

$$B\_{1,0}e^{i\vec{k}\_{m,n}\overline{x}} = -\frac{\text{i}aK\_Re^{i\mathcal{O}\_R}}{\sigma} \frac{1}{\cosh k\_{1,0}h\_1},\tag{31}$$

$$A\_{M,0}e^{-\mathrm{i}\vec{k}\_{M,0}\overline{\mathbf{x}}\_{M-1}} = -\frac{\mathrm{i}aK\_Te^{\mathrm{i}\partial\_T}}{\sigma}\frac{1}{\cosh k\_{M,0}\mathrm{l}\mathbf{t}\_M},\tag{32}$$

$$A\_{1,n} = 0 \text{ for } n = 1, 2, \dots, N\_\prime \tag{33}$$

and

$$B\_{M,n} = 0 \text{ for } n = 0, 1, \dots, N. \tag{34}$$

The other coefficients *Am*,*<sup>n</sup>* and *Bm*,*<sup>n</sup>* in Equation (20) can be determined using the matching conditions, Equations (8)–(11), at two adjacent shelves. The conservation of momentum, stemming from Equation (8), gives

$$
\left< \zeta\_{m,l}^{\text{inner}} \middle| \phi\_m \right> \Big|\_{x=x\_m} = \left< \zeta\_{m+1,l}^{\text{inner}} \middle| \phi\_{m+1} \right> \Big|\_{x=x\_m} \tag{35}
$$

where ζinner *<sup>m</sup>*,*l* (*z*) is the inner depth eigenfunction corresponding to *d*max *m* and *h*min *m* . For clarity, ζinner *<sup>m</sup>*,*l* (*z*) is defined by Equation (21) with the wavenumbers *km*,*<sup>n</sup>* obtained from Equations (16)–(19) as *dm* and *hm* replaced by *d*max *m* and *h*min *m* , respectively. Similarly, the conservation of mass, comes from Equations (9)–(11), yields the following equation

$$\left. \left\langle \frac{\partial \phi\_m}{\partial \mathbf{x}} \middle| \zeta\_{m,I}^{\text{outer}} \right\rangle \right|\_{\mathbf{x} = \mathbf{x}\_m} = \left\langle \frac{\partial \phi\_{m+1}}{\partial \mathbf{x}} \middle| \zeta\_{m,I}^{\text{outer}} \right\rangle \Big|\_{\mathbf{x} = \mathbf{x}\_m} \tag{36}$$

where ζouter *<sup>m</sup>*,*l* (*z*) is the outer depth function corresponding to *d*min *m* and *h*max *m* . In Equations (35) and (36), the subscripted indices go as *l* = 0, 1, ... , *N* and *m* = 1, 2, ... , *M* − 1. Additionally, it should be noted that Equations (35) and (36) are valid for all eight cases, as shown in Figure 2.

By using Equation (20), Equations (35) and (36) can be rewritten in the following forms

$$\begin{aligned} &\sum\_{m=0}^{N} \left( A\_{m,n} \underline{\boldsymbol{\xi}}\_{m,n}^{(1)}(\mathbf{x}\_{m}) + B\_{m,n} \underline{\boldsymbol{\xi}}\_{m,n}^{(2)}(\mathbf{x}\_{m}) \right) \underline{\boldsymbol{\zeta}}\_{m,l}^{\text{inner}}(\boldsymbol{\zeta}\_{m,n}) \\ &= \sum\_{m=0}^{N} \left( A\_{m+1,n} \underline{\boldsymbol{\xi}}\_{m+1,n}^{(1)}(\mathbf{x}\_{m}) + B\_{m+1,n} \underline{\boldsymbol{\xi}}\_{m+1,n}^{(2)}(\mathbf{x}\_{m}) \right) \underline{\boldsymbol{\zeta}}\_{m,l}^{\text{inner}}(\boldsymbol{\zeta}\_{m+1,n}) \end{aligned} \tag{37}$$

and

$$\begin{array}{l} \sum\_{n=0}^{N} \left( A\_{m,n} \frac{d\boldsymbol{\xi}\_{m,n}^{(1)}}{dx} (\boldsymbol{\chi}\_{m}) + B\_{m,n} \frac{d\boldsymbol{\xi}\_{m,n}^{(2)}}{dx} (\boldsymbol{\chi}\_{m}) \right) \Big\langle \zeta\_{m,n} \Big| \boldsymbol{\zeta}\_{m,l}^{\text{outer}} \right\rangle \\\ = \sum\_{n=0}^{N} \left( A\_{m+1,n} \frac{d\boldsymbol{\xi}\_{m+1,n}^{(1)}}{dx} (\boldsymbol{\chi}\_{m}) + B\_{m+1,n} \frac{d\boldsymbol{\xi}\_{m+1,n}^{(2)}}{dx} (\boldsymbol{\chi}\_{m}) \right) \Big\langle \zeta\_{m+1,n} \Big| \boldsymbol{\zeta}\_{m,l}^{\text{outer}} \Big\rangle. \end{array} \tag{38}$$

Subsequently, it can be found that Equations (30), (33), (34), (37), and (38) are 2*M*(*N*+1) linear equations, that can be used to solve the 2*M*(*N*+1) unknowns, *Am,n* and *Bm,n*. Furthermore, Equations (37) and (38) can be reduced to the original equations of EMM for normal water-wave-scattering without structures where α = 0 and *dm* = 0 [28,35]. In this study, the SuperLU is used to solve the resulting sparse system of linear equations [52]. After the unknowns, *Am,n* and *Bm,n*, are solved, the reflection and transmission coefficients can be obtained by Equations (31) and (32), respectively. This completes the solution procedure for the EMM.

**Figure 2.** Schematics for eight different situations of shelves separated by abrupt connections.
