## 2 Related Works
### The Updated Algorithm Considering Vertical Turbulent Wind Field
\begin{array}{cccc}1&e^{i\theta_{12}(\omega) }&\cdots&e^{i\theta_{1n}(\omega)}\\ e^{i\theta_{21}(\omega)}&1&\cdots&e^{i\theta_{2n}(\omega)}\\ \vdots&\vdots&\ddots&\vdots\\ e^{i\theta_{n1}(\omega)}&e^{i\theta_{n2}(\omega)}&\cdots&1\end{array} \right] \tag{28}\]  
where \(\otimes\) denotes the element-wise multiplication operator between two matrices. Similarly, \(\Gamma(\omega)\) can be given by:  
\[\mathbf{\Gamma}(\omega)=\big{|}\mathbf{\Gamma}(\omega)\big{|}\otimes\mathbf{ \Theta}(\omega) \tag{29}\]  
where \(\big{|}\Gamma(\omega)\big{|}\) represent the lagged coherence matrix [33]. According to Equation (8), \(\big{|}\mathbf{S}(\omega,t)\big{|}\) can be factorized as:  
\[\big{|}\mathbf{S}(\omega,t)\big{|}=\mathbf{\Lambda}(\omega,t)\mathbf{\Lambda }^{T}(\omega,t) \tag{30}\]  
where \(\Lambda(\omega,t)\) is lower triangular matrix. The following relations have been proved when the phases admit the requirement that \(\theta_{jm}(\omega)=\theta_{jk}(\omega)-\theta_{km}(\omega)\)[32]:  
\[H_{jj}(\omega,t)=\Lambda_{jj}(\omega,t)\text{; }H_{jm}(\omega,t)=\Lambda_{jm}( \omega,t)e^{i\theta_{jm}(\omega)}\text{ }(j>m) \tag{31}\]  
\[\theta_{jm}(\omega)=\vartheta_{jm}(\omega) \tag{32}\]Similar to Equation (23), \(\left|\mathbf{S}(\omega,t)\right|\) can be rewritten as:  
\[\left|\mathbf{S}(\omega,t)\right|=\mathbf{D}(\omega,t)\mathbf{B}(\omega)\mathbf{ B}^{T}(\omega)\mathbf{D}^{T}(\omega,t) \tag{33}\]  
where \(\mathbf{B}(\omega)\) is the lower triangular matrix derived from \(\left|\Gamma(\omega)\right|\) through Cholesky decomposition. That is:  
\[\left|\mathbf{\Gamma}(\omega)\right|=\mathbf{B}(\omega)\mathbf{B}^{T}(\omega) \tag{34}\]  
where  
\[\mathbf{B}(\omega)=\left[ \begin{array}{cccc}\beta_{11}(\omega)&0&\cdots&0\\ \beta_{21}(\omega)&\beta_{22}(\omega)&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ \beta_{n1}(\omega)&\beta_{n2}(\omega)&\cdots&\beta_{nn}(\omega)\end{array} \right] \tag{35}\]  
According to Equations (30) and (33), the following relationship can be obtained:  
\[\mathbf{\Lambda}(\omega,t)=\mathbf{D}(\omega,t)\mathbf{B}(\omega) \tag{36}\]  
Based on Equations (31), (32) and (36), the components of the \(\mathbf{x}(t)\) can be generated using:  
\[x_{j}(t)=\text{Re}\Bigg{(}2\sum_{m=1}^{j}\sum_{l=1}^{N}\sqrt{S_{jj}\big{(} \omega_{l^{\prime}}t\big{)}\Delta\omega}\beta_{jm}\Big{(}\omega_{l}\Big{)}e^{ i(\omega_{l}t-\theta_{jm}(\omega_{l})+\Phi_{ml})}\Bigg{\}} \tag{37}\]  
Taking advantage of the NMF, the square root of \(S_{ii}(\omega,t)\) can be factorized and given as:  
\[\sqrt{S_{jj}(\omega,t)}\approx\sum_{q=1}^{P}w_{q}^{jj}(t)v_{q}^{jj}(\omega) \tag{38}\]  
Substituting Equation (38) into Equation (37), the updated NMF-based non-stationary SRM can be rewritten as:  
\[x_{j}(t)=\text{Re}\Bigg{\{}2\sum_{m=1}^{j}\sum_{q=1}^{P}w_{q}^{jj}(t)\sum_{l=1 }^{N}v_{q}^{jj}\big{(}\omega_{l}\big{)}\beta_{jm}\Big{(}\omega_{l}\Big{)} \sqrt{\Delta\omega}e^{i(\omega_{l}t-\theta_{jm}(\omega_{l})+\Phi_{ml})} \Bigg{\}} \tag{39}\]  
Or  
\[x_{j}(t)=\text{Re}\Bigg{\{}2\sqrt{\Delta\omega}\sum_{q=1}^{P}w_{q}^{jj}(t) \sum_{l=1}^{N}v_{q}^{jj}\big{(}\omega_{l}\big{)}e^{i\omega_{l}t}\sum_{m=1}^{j }\beta_{jm}\big{(}\omega_{l}\big{)}e^{i(\Phi_{ml}-\theta_{jm}(\omega_{l}))} \Bigg{\}} \tag{40}\]  
It is clear that the FFT technique can be introduced for the second summation in Equation (40), which can enhance the simulation efficiency significantly. In addition, various auto-EPSDs and the complex coherence functions can be considered in Equation (40) for the simulation of both horizontal and vertical non-stationary turbulent wind fields for engineering structures.