## 3 Optimization of the Phase Screen Model  
Modeling errors will be induced into the phase screens because of the numerical sampling. We propose using the subharmonic method [28] and the strict sampling constraint method to compensate for the insufficient discrete sampling. First, the sampling points of low spatial frequencies are increased by adding a subharmonic phase screen to the original screen. A subharmonic phase screen is a sum of \(N_{p}\)-layers of screens given by:  
\[\varphi_{LF}(x,y)=\sum_{p=1}^{N_{p}}\sum_{n=-1}^{1}\sum_{m=-1}^{1}c_{n,m}\exp \left[j\left(k_{xn}x+k_{ym}y\right)\right] \tag{10}\]  
where each value of \(p\) corresponds to a different \(3\times 3\) grid of frequency, and the frequency grid spacing is \(\Delta f_{p}=1/(3^{p}\ L)\). Moreover, the geometry (taking the \(y\)-\(z\) plane for example) of the diffraction process is illustrated in Figure 2, where \(D_{1}\) and \(D_{2}\) denote the interested observation region on the screens, and \(\alpha\) is the deviation angle of the propagation beam; the maximum angle is reckoned as:  
\[\sin\alpha_{max}=\frac{D_{1}+D_{2}}{2\Delta z} \tag{11}\]According to the Nyquist sampling theorem, the grid spacing \(d\) of the numerical simulation model should be constrained as \(d<1/(2f_{max})\), where \(f_{max}\) is the maximum spatial frequency of interest and satisfies \(f_{max}=\sin\alpha_{max}/\lambda\). Therefore, the upper bound of the grid spacing \(d\) can be derived as:  
\[d_{max}=\frac{\lambda\Delta\mathrm{z}}{D_{1}+D_{2}} \tag{12}\]  
According to Equation (12), strict constraint on \(d\) can be applied to improve the accuracy of the model since the larger the observation region (\(D_{1},D_{2}\)), the more complete the information. The joint application of the subharmonic method and the strict sampling constraint method provides flexibility for the model optimization, consistent with the analytical theory. The phase structure function (PSF) is employed as a statistical reference to verify the accuracy of the established phase screen model. The theoretical expression of the underwater turbulence PSF is [29]:  
\[D_{\varphi}^{theory}(r)=3.603\times 10^{-7}k^{2}z\varepsilon^{-1/3}\chi_{T}r^{ 5/3}\Big{(}1.16-2.235/\omega+1.119/\omega^{2}\Big{)} \tag{13}\]  
where \(r\) is the separation distance between two points on the phase front. Meanwhile, PSF generated from the existing phase screens is defined as [30]:  
\[D_{\varphi}(r)=\langle\big{[}\varphi\big{(}r^{\prime}+r\big{)}-\varphi\big{(} r^{\prime}\big{)}\big{]}^{2}\rangle \tag{14}\]  
where \(r^{\prime}\) represents the coordinate on a phase screen. A statistical average of 1000 temporal-coherent phase screens is carried out and then compared with the theoretical PSF. The channel parameters are set as follows: \(\lambda=450\,\mathrm{nm}\), \(\Delta\mathrm{z}=20\,\mathrm{m}\), \(\varepsilon=10^{-4}\,\mathrm{m}^{2}/\mathrm{s}^{3}\), \(\chi_{\mathrm{T}}=10^{-6}\,\mathrm{K}^{2}/\mathrm{s}\), \(\omega=-2\), \(\eta=1\,\mathrm{mm}\), \(a=0.0405\,\mathrm{m}^{-1}\), and \(b=0.0025\,\mathrm{m}^{-1}\). We also define the matching degree between the model and theory by normalized 2-norm as:  
\[M=1-\frac{\|D_{\varphi}^{theory}-D_{\varphi i}\|_{2}}{\|D_{\varphi}^{theory}- D_{\varphi 0}\|_{2}} \tag{15}\]  
where \(D_{\varphi 0}\) is the initial PSF of the model without optimization and \(D_{\varphi i}\) is each PSF obtained by adjusting the optimization conditions. The range of \(M\) is (\(-\infty\), 1); the closer the value is to 1, the better the matching degree. Figure 3 shows the optimization performance of the subharmonic method. The PSF curve becomes closer to the dotted line representing the theoretical value, with the subharmonic layer number \(N_{p}\) increasing from 0 to 5. However, at \(N_{p}=6\), the low-frequency sampling is overcompensated, resulting in the deviation from the theoretical curve again. The PSF matching degrees with \(N_{p}=0\), 3, 5, and 6 are 0, 1.07%, 48.09%, and \(-34.68\%\), respectively.