## 2 Methodology of the Forward Modelling  
### Approximation of AVA Seismic Data Simulation  
The propagation matrix method requires a constant slowness \(s_{1}\), as indicated by the equation \(s_{1}=\frac{\sin(\theta)}{V_{\rm P}}\), where the P-wave velocity variations \(V_{\rm P}\) in the subsurface medium cause changes in the incident angle of the seismic waves. \(\theta\) represents the incident angle of the P-wave seismic waves. In this study, the propagation matrix theory was modified to simulate AVA pre-stack seismic responses more accurately. It is assumed that the seismic wave amplitudes are constant within the geological layers, and the reflection interfaces exhibit continuous stress-strain and a constant incident angle for the P-wave. Specifically, if there are \(nz\)-1 reflection interfaces in the subsurface, the incident wave at the \(iz\)-th interface with the incidence angle \(\theta\) equals the transmitted wave at the same angle at the (\(iz\)\(-\)1)-th interface. These modifications result in a new propagation matrix theory for simulating AVA pre-stack seismic responses. According to [PERSON] [27,30], the relationship between the particle velocity/stress array and the amplitude of seismic waves can be expressed by the following equation:  
\[ \begin{pmatrix}v_{1}\\ v_{3}\\ \sigma_{33}\\ \sigma_{13}\end{pmatrix} =\mathbf{T}(z) \begin{pmatrix}U_{\rm P}^{-}\\ U_{\rm S}^{-}\\ U_{\rm P}^{+}\\ U_{\rm S}^{+}\end{pmatrix} . \tag{1}\]  
where \(v\) denotes the vibration velocity of the particle; \(\sigma\) represents the stress component of the particle; subscripts 1 and 3 represent the x-direction and z-direction, respectively; \(U^{-}\) and \(U^{+}\) represent the amplitude of the upward and downward wave, respectively; P and S represent the P-wave and S-wave, respectively. \(\mathbf{T}\) can be expressed as follows:  
\[\mathbf{T}(z)=\mathbf{F}_{iz,1}\mathbf{E}_{iz,1}. \tag{2}\]  
In this context, \(iz\) denotes the index of the \(iz\)-th reflector interface. Furthermore, \(z\) represents the depth of the calculation point relative to this specific reflector interface within the geological structure. Matrices \(\mathbf{F}_{iz,1}\) and \(\mathbf{E}_{iz,1}\) can be expressed in the following forms:  
\[\mathbf{F}_{iz,1}={\rm i}\omega \begin{pmatrix}\beta_{\rm P}^{iz,1}&\beta_{\rm S}^{ iz,1}&\beta_{\rm P}^{iz,1}&\beta_{\rm S}^{iz,1}\\ -\frac{\beta_{\rm P}^{iz,1}}{-\beta_{\rm S}^{iz,1}}&\beta_{\rm P}^{iz,1}&\beta _{\rm S}^{iz,1}\\ -Z_{\rm P}^{iz,1}&-Z_{\rm P}^{iz,1}&-Z_{\rm P}^{iz,1}&-Z_{\rm P}^{iz,1}\\ W_{\rm P}^{iz,1}&W_{\rm S}^{iz,1}&-W_{\rm P}^{iz,1}&-W_{\rm S}^{iz,1}\end{pmatrix} , \tag{3}\]  
\[\mathbf{E}_{iz,1}={\rm diag}\Big{[}{\rm e}^{i\omega_{\rm S}^{iz,1}z}\cdot{\rm e} ^{i\omega_{\rm S}^{iz,1}z}\cdot{\rm e}^{-i\omega_{\rm S}^{iz,1}z}\cdot{\rm e} ^{-i\omega_{\rm S}^{iz,1}z}\cdot{\rm e}^{-i\omega_{\rm S}^{iz,1}z}\Big{]}. \tag{4}\]  
In this context, \(\omega\) represents the angular frequency, where \(\omega=2\pi f\), with \(f\) denoting the computation frequency. Furthermore, i is used as the symbol for the imaginary unit. The calculations of \(\beta\), \(\xi\), \(W\), and \(Z\) will be introduced in Equations (9)-(11) below. The amplitude relationship on both sides of the reflection interface can be expressed using Equation (1) as follows:  
\[\Big{(}\mathbf{F}_{iz,1}\mathbf{E}_{iz,1}\bigg{(}\frac{1}{2}z\bigg{)}\Big{)}\Big{\backslash} \Big{(}\mathbf{F}_{iz+1,2}/\mathbf{E}_{iz+1,2}\Big{(}\frac{1}{2}z\Big{)}\Big{)} \begin{pmatrix}U_{iz+1,\rm P}^{-}\\ U_{iz+1,\rm S}^{-}\\ U_{iz+1,\rm P}^{+}\\ U_{iz+1,\rm P}^{+}\\ U_{iz+1,\rm S}^{+}\end{pmatrix} = \begin{pmatrix}U_{iz,\rm P}^{-}\\ U_{iz,\rm P}^{-}\\ U_{iz,\rm P}^{-}\\ U_{iz,\rm P}^{+}\\ U_{iz,\rm S}^{+}\end{pmatrix} .