*Appendix B.1. Horizontal Layer Model*

We start with a simple layered model with *N* layers and *N* − 1 interfaces. The P-wave velocity map is used to describe the model. For each column in the area, the points are filled with preset mediums using:

$$M(X,Y) = m\_{\text{ul}} \; , \; y\_n < Y \le y\_{n+1} \; , \; n = 0, 1, 2 \dots, N \tag{A19}$$

where *M*(*X*,*Y*) is the elastic parameter of the point (*X*,*Y*), *mk* is the preset parameter of medium *u*, *yn* and *yn*+<sup>1</sup> are respectively the depth of the *n*th and (*n* + 1)th interface. The three-layer (two-interface) model is shown in Figure A3.

**Figure A3.** Horizontal layered model.

*Appendix B.2. Fold*

A Gaussian curve is used to form folds in the model because a Gaussian curve can satisfactorily present folds when the PML layer and fault dislocation are cut (Figure A4).

$$y\_{\mathcal{S}}(x) = a \mathbf{e}^{-(x-b)^2/2c^2} \tag{A20}$$

where *a* is the fold depth, *b* is the horizontal location of the fold and *c* is a parameter controlling the width of the fold.

**Figure A4.** Gaussian curve used to control the fold shape.

The interfaces between the layers are often not totally smooth. The curve *yg*(*x*) is used as a baseline to create fluctuations. The fluctuation *yf* is composed of random floating values with 0 as the mean. Now, the interfaces *yi* between layers (Figure A5) can be obtained by Equation (A20). Since the Gaussian curve and fluctuations are randomly generated, the interfaces can be either parallel or nonparallel (Figure A5).

$$y\_i(\mathbf{x}) = y\_n + y\_\mathcal{J}(\mathbf{x}) + y\_f(\mathbf{x})\tag{A21}$$

**Figure A5.** Interface with fluctuations.

**Figure A6.** Fold model (**a**) without and (**b**) with random fluctuations.

### *Appendix B.3. Fault*

The fault is determined by the fault line and dislocation distance. The fault line controls the position and inclination angle of the fault. The dislocation distance controls the amount of fault displacement. We select a point (*x*0, *y*0) in the area as the reference point and use a point-slope linear equation with a slope of *k*.

$$\frac{\mathcal{Y} - \mathcal{y}\_0}{\mathcal{X} - \mathcal{x}\_0} = k \tag{A22}$$

When the horizontal displacement Δ*x* is determined, the fault can be generated as:

$$M(X,Y) = M(\mathbf{X} + \Delta \mathbf{x}, Y + k\Delta \mathbf{x}), \ X < (Y - y\_0)/k + \ x\_{0\prime}(X, Y) \in \Omega \tag{A23}$$

Equation (A23) shows that the points on the left side of the fault line move along the fault line by Δ*x* in the horizontal direction and *k*Δ*x* in the vertical direction. Fault generation is also a random procedure. To ensure the stability of random processing and the existence of sufficient space for the PML boundary after the fault slips, a 10 m-thick layer is added around the model. As shown in Figure A7, blue represents the model area, orange shows the PML boundary, red indicates the layer reserved for fault dislocation, and light blue shows the simulation area after fault generation.

**Figure A7.** Fault generation illustration. (**a**) Model before fault displacement with the fault line and reference point. (**b**) Generated fault model.

### *Appendix B.4. Cave*

Caves filled with low-velocity materials (*mc*) are added into the model as anomalies by simply controlling the cave center location (*xc*, *yc*) and radius *rc*:

$$M(\mathbf{X}, \mathbf{Y}) = m\_{\varepsilon \prime} \left( X - \mathbf{x}\_{\varepsilon} \right)^{2} + \left( \mathbf{Y} - y\_{\varepsilon} \right)^{2} \le r\_{\varepsilon}^{2} \text{ } (\mathbf{X}, \mathbf{Y}) \in \Omega \tag{A24}$$

Now, we can obtain a model with horizontal parallel layers, folds, random fluctuations, faults, and caves incorporated (Figure A8).

**Figure A8.** Model with folds, random fluctuations, faults, and caves.
