Method and system for simulating the interaction of seafloor topography with ocean currents

Abstract

The application provides a simulation method and system for the interaction between seabed topography and ocean current, comprising: obtaining high-resolution seabed topography data, generating an adaptive grid structure, and constructing fine topography description data; generating interaction data between the ocean current field and the topography by using a turbulence closure model according to the fine topography description data; simulating a nonlinear process to generate nonlinear process characteristic data according to the interaction data; generating ocean mixing process data by using a particle tracking method according to the nonlinear process characteristic data; optimizing turbulence parameters to generate optimized turbulence mixing data according to the ocean mixing process data; and updating the simulation of the ocean current field to generate high-precision ocean current field distribution data according to the optimized turbulence mixing data. The application can improve the simulation accuracy of ocean dynamic processes under complex topography conditions and provide important support for ocean environment prediction and resource development.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of information technology, and in particular to a simulation method for the interaction between seabed topography and ocean currents. BACKGROUND

[0002] Ocean numerical simulation is a core tool for studying ocean dynamic processes and predicting ocean environmental changes, which is of great significance for climate change, ocean ecological protection and disaster warning. Through the construction of mathematical models, scientists can simulate the spatiotemporal evolution of ocean currents, temperature, salinity and other elements. However, existing methods have significant limitations in dealing with complex ocean environments. Traditional models often rely on simplified parameterization schemes, resulting in insufficient accuracy in simulating the interaction of seabed topography, coastlines and water masses, especially in expressing strong nonlinear processes such as density fronts or internal wave breaking, which makes it difficult to capture the details of real physical processes. In addition, the simulation of momentum and energy exchange at the sea-air and sea-ice interfaces often ignores the complexity of multi-medium coupling, which weakens the predictive ability of the model for real ocean systems.

[0003] These limitations reveal the core challenges in the field of ocean numerical simulation. First, the simulation of complex topography and ocean current interaction requires high-precision boundary layer parameterization schemes, but current methods struggle to balance computational efficiency and physical realism. This leads to models often failing to accurately reflect the impact of topography on the flow field when simulating the seabed boundary layer, which in turn weakens the accurate depiction of water mass movement. The resulting technical problem is how to achieve fine simulation of strong nonlinear processes such as density fronts or internal wave breaking under limited computing resources, which are crucial for ocean mixing and energy transfer. Further, accurately simulating these nonlinear processes requires improved coupling mechanisms at the sea-air and sea-ice interfaces, as the momentum and energy exchange at the interface directly affects the dynamic behavior of water masses.

[0004] Therefore, how to develop high-precision seabed boundary layer parameterization schemes, improve numerical representation of water mass collision mixing, and optimize the simulation of momentum and energy exchange at the sea-air and sea-ice interfaces, has become a key issue for improving the predictive ability of ocean numerical simulation. SUMMARY

[0005] The present application provides a simulation method for the interaction between seabed topography and ocean currents, mainly comprising:

[0006] The method comprises the following steps: acquiring high-resolution seafloor topographic data, generating an adaptive grid structure, constructing fine topographic description data, generating interaction data of the seafloor topographic data and a sea current field by using a turbulence closure model according to the fine topographic description data, simulating a nonlinear process according to the interaction data to generate nonlinear process characteristic data, generating marine mixing process data by using a particle tracking method according to the nonlinear process characteristic data, optimizing turbulence parameters according to the marine mixing process data to generate optimized turbulence mixing data, and updating the simulation of the sea current field according to the optimized turbulence mixing data to generate high-precision sea current field distribution data.

[0007] Further, the method of acquiring high-resolution seafloor topographic data, generating an adaptive grid structure, and constructing fine topographic description data comprises the following steps: acquiring a three-dimensional coordinate point set from a high-resolution seafloor topographic data set; calculating a topographic gradient according to the three-dimensional coordinate point set; if the topographic gradient is greater than a preset threshold, using a quadtree algorithm to encrypt a grid; and generating fine topographic description data by using a bilinear interpolation method through the encrypted grid.

[0008] Further, the method of generating interaction data of the seafloor topographic data and a sea current field by using a turbulence closure model according to the fine topographic description data comprises the following steps: extracting a topographic feature point set from the fine topographic description data; calculating a turbulence mixing coefficient by using a turbulence closure model; extracting a sea current field shear feature from the turbulence mixing coefficient; and generating the interaction data of the sea current field and the topography by using a principal component analysis method to reduce the dimension of the interaction feature vector.

[0009] Further, the method of simulating a nonlinear process according to the interaction data to generate nonlinear process characteristic data comprises the following steps: acquiring a nonlinear feature distribution from the interaction data; calculating a nonlinear term energy by using a high-order spectrum method; if the nonlinear term energy is greater than a preset threshold, introducing a sub-grid model; and generating the nonlinear process characteristic data by using a finite volume method through the sub-grid model.

[0010] Further, the method of generating marine mixing process data by using a particle tracking method according to the nonlinear process characteristic data comprises the following steps: extracting a velocity field data from the nonlinear process characteristic data; calculating a particle trajectory distribution by using a Lagrangian particle tracking method; calculating a velocity gradient field according to the particle trajectory distribution; and extracting a density front and internal wave breaking feature from the velocity gradient field to generate the marine mixing process data.

[0011] Further, the optimizing the turbulent flow parameter according to the ocean mixing process data to generate optimized turbulent mixing data comprises: acquiring a turbulent flow feature from the ocean mixing process data; calculating an initial turbulent viscosity coefficient by using a turbulent closure model; if the initial turbulent viscosity coefficient deviates from the turbulent flow feature by more than a preset threshold, updating the turbulent viscosity coefficient by using an iterative optimization method to generate the optimized turbulent mixing data.

[0012] The technical scheme provided by the embodiment of the application can have the following beneficial effects:

[0013] The application discloses a fine simulation method for ocean dynamic processes based on high-resolution seafloor topographic data. The method adopts an adaptive grid division technique, dynamically adjusts the grid density according to the topographic gradient, and constructs a seafloor boundary layer parameterization scheme. The terrain-induced turbulent mixing is calculated by using a Mellor-Yamada turbulent closure model, and a high-order spectral method is used to simulate nonlinear processes. The Lagrangian particle tracking method is used to extract the spatiotemporal distribution of density fronts and internal wave breaking, and the turbulent parameterization scheme is optimized. Finally, the application applies a three-dimensional non-hydrostatic ocean model to update the simulation of the ocean current field, evaluates the model accuracy through cross-validation, and feeds back the optimization of grid division and the turbulent closure model. The method significantly improves the simulation accuracy of ocean dynamic processes under complex topographic conditions, and provides important support for ocean environment prediction and resource development. BRIEF DESCRIPTION OF DRAWINGS

[0014] Figure 1 The flowchart of the simulation method for the interaction between seafloor topography and ocean currents.

[0015] Figure 2 The flowchart of generating adaptive grids. DETAILED DESCRIPTION

[0016] In order to better understand the content of the application, the application will be described in detail in combination with the drawings and embodiments. The specific embodiments described herein are only used to explain the related application, and are not limited to the application. In addition, it should be noted that, in order to facilitate description, only the parts related to the application are shown in the drawings.

[0017] As Figure 1 The simulation method for the interaction between seafloor topography and ocean currents can specifically comprise:

[0018] Step S101, please refer to Figure 2Data is obtained from a high-resolution seafloor topography dataset, and an adaptive meshing technique is used to adjust the grid density according to the terrain gradient. For example, a gradient threshold of 0.01 is set, with a unit of height change per kilometer. If the gradient is greater than the threshold, the grid is encrypted to the minimum unit of 0.1 kilometers. At this time, fine terrain description data is generated.

[0019] The original terrain data obtained from the high-resolution seafloor topography dataset is in the format of a digital elevation model. Through raster data analysis, a set of three-dimensional coordinate points is obtained. According to the set of three-dimensional coordinate points, the height change per kilometer is calculated. The gradient calculation formula can be used where H1 and H2 are the heights of adjacent points, and D is the distance. The terrain gradient distribution is obtained. If the terrain gradient is greater than the preset gradient threshold, the quadtree algorithm is used for grid encryption in the corresponding area to generate an adaptive grid structure. Through the adaptive grid structure, the terrain height in each grid cell is calculated using a bilinear interpolation method to generate fine terrain description data. Specifically, first, the original elevation data is extracted from the high-resolution seafloor topography dataset, with a data accuracy of 0.5 kilometers per pixel. The data is resampled to a resolution of 0.2 kilometers using a bicubic interpolation algorithm to ensure the integrity of the terrain details. Based on the resampled data, the terrain gradient value of each grid cell is calculated, and the Sobel operator is used for gradient detection. The horizontal direction kernel is [-1 0 1; -2 0 2; -1 0 1], and the vertical direction kernel is its transpose. The gradient amplitude of each point is obtained through convolution operation. The gradient threshold is set to 0.01 (i.e., 1 meter of height change per kilometer). When a region with a gradient value exceeding the threshold is detected, the quadtree subdivision algorithm is started to recursively subdivide the region. Each subdivision divides the current grid into 4 subgrids until the minimum unit size of 0.1 kilometers or the gradient value is below the threshold.

[0020] During the grid encryption process, the newly generated grid nodes are processed using the constrained Delaunay triangulation algorithm to ensure the topological consistency between adjacent grids. The Kriging interpolation method (semi-variogram function using Gaussian model, with a range of 0.3 kilometers) is used to assign elevation values to the new nodes. The final fine terrain data uses a hierarchical storage structure, with the base layer retaining 0.2 kilometer grids and the detail layer storing 0.1 kilometer encrypted grids. Spatial indexing (R-tree structure with a node capacity of 50) is used to achieve fast retrieval and visualization rendering of multi-resolution data. For gradient mutation regions (such as trench edges), the system automatically marks them as feature boundaries and prioritizes super-resolution processing at 0.05 kilometers in subsequent analysis. Convolutional neural networks (U-Net architecture with 5 layers of downsampling and 5 layers of upsampling) are used to enhance and reconstruct the local terrain.

[0021] The application realizes local grid refinement through a quadtree encryption algorithm triggered by a terrain gradient threshold (such as 0.01 / km), generates 0.1km or even 0.05km high-resolution grids in areas with sudden changes in seabed topography (such as trenches and seabed mountains), and solves the problem of poor adaptability of traditional uniform grids to complex topography. By bilinear interpolation in the encrypted grid, the original elevation details are preserved (such as complete reconstruction of 0.5km resolution data), and the computational complexity caused by global encryption is avoided, optimizing the computational efficiency. For areas with sudden changes in gradient (such as the edges of trenches), super-resolution processing (0.05km grid) is used in combination with a CNN model to reconstruct the local terrain, improving the physical realism of the terrain description and providing high-precision boundary conditions for subsequent turbulence calculations, enhancing the topographic features.

[0022] Step S102, according to the refined terrain description data obtained in step S101, construct a seabed boundary layer parameterization scheme, use the Mellor-Yamada turbulence closure model to calculate the terrain-induced turbulent mixing coefficient, and generate seabed current field interaction expression data.

[0023] From the refined terrain description data, extract the terrain feature point set, use the gradient-based feature detection method to calculate the terrain surface height change rate, and obtain the terrain feature distribution. According to the terrain feature distribution, use the Mellor-Yamada turbulence closure model to calculate the turbulent mixing coefficient, and through the formula Km = l^2Sf(Ri), where Km is the turbulent mixing coefficient, l is the mixing length, S is the flow field shear rate, Ri is the Richardson number, and f(Ri) is the stability function, to obtain the turbulent mixing coefficient distribution. From the turbulent mixing coefficient distribution, extract the seabed current field shear feature, use the finite difference method to calculate the seabed current velocity gradient, and through the formula where ui is the seabed current velocity, z is the vertical coordinate, and Δz is the grid spacing, to obtain the seabed current field shear distribution. If there is an area in the seabed current field shear distribution with a shear rate exceeding a preset threshold, then for the area, use the k-ε turbulence model to calculate the local turbulent kinetic energy, through the formula k = 0.5(u'^2 + v'^2 + w'^2), where k is the turbulent kinetic energy, u', v', and w' are the turbulent velocity components, to obtain the local turbulent kinetic energy distribution. According to the local turbulent kinetic energy distribution, calculate the terrain-induced turbulent enhancement effect, use the bilinear interpolation method to interpolate the turbulent kinetic energy to the terrain feature point set, and obtain the terrain interaction enhancement data. From the terrain interaction enhancement data, extract the seabed current field and terrain interaction features, use the principal component analysis method to reduce the dimension of the interaction feature vector, and obtain the interaction expression feature set. According to the interaction expression feature set, generate structured data of the seabed current field and terrain interaction, use the grid mapping method to map the feature set to a three-dimensional coordinate grid, and obtain the seabed current field interaction distribution data.

[0024] In the construction of the seabed boundary layer parameterization scheme based on refined topographic description data, the Mellor-Yamada turbulence closure model is adopted to generate the expression data of the interaction between the sea current field and the topography by calculating the topography-induced turbulent mixing coefficient. First, the elevation information is extracted from the refined topographic data, and the topographic curvature is calculated using the finite difference method, with the curvature threshold set to 0.005 per square kilometer. When the curvature exceeds the threshold, it is marked as a significant area of turbulent mixing. In the Mellor-Yamada model, the topographic curvature is introduced as a correction factor for the turbulent mixing coefficient. The second-order moment closure method is used to calculate the turbulent kinetic energy (TKE) and the turbulent dissipation rate (ε), with the initial value of TKE set to 0.01 square meters per square second and the initial value of ε set to 0.001 square meters per cubic second. By iteratively solving the turbulence closure equation set, the spatial distribution of the turbulent mixing coefficient is obtained, with the mixing coefficient range controlled between 0.1 and 10 square meters per second. Further, the turbulent mixing coefficient is coupled with the sea current field, and the Navier-Stokes equation is used to describe the sea current motion, with a horizontal resolution of 0.2 kilometers and a vertical stratification of 10 layers, each with a thickness of 50 meters. Through numerical simulation, the interaction between the sea current field and the topography is calculated to generate the spatial distribution data of the flow velocity, flow direction, and turbulent intensity. For complex topographic areas such as trenches and submarine mountains, the adaptive grid refinement technique is used to increase the grid resolution to 0.05 kilometers to ensure the accuracy of the simulation results. Finally, the simulation data is resampled to a uniform grid using a spatial interpolation algorithm (such as the inverse distance weighted method with a weight exponent of 2) to generate comprehensive expression data of the interaction between the sea current field and the topography, which can be used for subsequent ocean dynamics research.

[0025] Based on the Mellor-Yamada turbulence closure model and the principal component analysis (PCA) processing flow, the topographic curvature is used to correct the turbulent mixing coefficient, solving the problem of ignoring topography-induced turbulence in traditional models and dynamically calibrating the turbulent mixing coefficient. PCA is used to reduce the dimensionality of the sea current shear characteristics (such as ), extract the dominant interaction mode (such as the coupling relationship between topographic slope and turbulent energy), and reduce the computational complexity of multi-physical field coupling. The Richardson number (Ri) dependent stability function f(Ri) is introduced to avoid the occurrence of unstable turbulent solutions and improve the robustness of the model in strong mixing areas.

[0026] Step S103, according to the expression data of the interaction between the sea current field and the topography, simulate the nonlinear process, use high-order spectral method to decompose the nonlinear term of density front and internal wave, if the energy of nonlinear term exceeds 0.001 square meter per square second, introduce sub-grid model, minimum unit 0.01 kilometers, generate nonlinear process detail feature data.

[0027] Nonlinear feature distribution is obtained from the interaction of the current field and the terrain expression data, high-order spectrum method is adopted, spectrum function is calculated through the formula F(k)=∫u(x)e^(-ikx)dx, wherein F(k) is the spectrum function, u(x) is the velocity field, k is the wave number, and x is the spatial coordinate, and the spectrum function distribution is obtained. According to the spectrum function distribution, the energy of the nonlinear term is calculated, and the energy distribution of the nonlinear term is obtained through the formula E=∫|F(k)|^2dk, wherein E is the energy, F(k) is the spectrum function, and k is the wave number. If there is a region exceeding the preset threshold in the energy distribution of the nonlinear term, a sub-grid model is introduced for the region, a finite volume method is adopted, and the discrete velocity field distribution is obtained through the formula Discretization, wherein u is the velocity, t is the time, p is the pressure, ρ is the density, and ν is the viscosity coefficient. The density front feature is extracted from the discrete velocity field distribution, the principal component analysis method is adopted, the velocity field matrix is decomposed through the formula X=UΣV^T, wherein X is the velocity field matrix, U is the left singular vector, Σ is the singular value, and V is the right singular vector, and the density front feature distribution is obtained. According to the density front feature distribution, the internal wave dynamic feature is calculated, and the internal wave dynamic feature distribution is obtained through the formula ω=√(N^2k_h^2 / (k_h^2+

[0028] k_z^2)), wherein ω is the internal wave frequency, N is the Brunt-Vaisala frequency, k_h is the horizontal wave number, and k_z is the vertical wave number. Nonlinear process detail features are extracted from the internal wave dynamic feature distribution, a bilinear interpolation method is adopted, and the features are interpolated to the grid points through the formula f(x,y)=(1-x)(1-y)f(0,0)+x(1-y)f(1,0)+(1-x)yf(0,1)+xyf(1,1), wherein f is the feature value, x and y are interpolation coordinates, and the nonlinear process detail feature data is obtained. According to the nonlinear process detail feature data, a structured feature set is generated, a grid mapping method is adopted, and the nonlinear process structured feature distribution is obtained through the formula G(i,j,k)=f(xi,yj,zk), wherein G is the grid feature, f is the feature value, xi,yj,zk are grid coordinates.

[0029] The present embodiment uses the interaction of the ocean current field and the terrain expression data as input, adopts the high-order spectrum method to decompose the nonlinear terms of the density front and internal wave, sets the Fourier mode truncation wave number to 256, the spatial resolution to 0.1 km, and the time step to 10 seconds. When calculating the nonlinear interaction term in the spectral space, if the energy spectrum density exceeds 0.001 square meter per square second, the sub-grid model processing mechanism is triggered, the dynamic adaptive threshold technology is adopted, the local grid is encrypted to 0.01 km, and the eddy viscosity coefficient is introduced for correction, with the value range of 0.05 to 0.5 square meters per second. For the energy concentrated nonlinear wave packet, the WKB approximation method is used to track the phase change, the phase resolution is set to 0.01 radian, and the weakly compressible fluid equation is coupled, with the Poisson ratio set to 0.3 and the sound speed gradient threshold set to 0.01 per second. In the sub-grid model, the high-frequency turbulent signal is extracted by using wavelet transform, the Daubechies-6 wavelet basis is selected, the decomposition scale is 5 layers, and the energy reconstruction error is controlled within 5%. The nonlinear coupled equation set is solved by the implicit iteration method, the relaxation factor is 0.7, the residual convergence standard is set to 1e-6, and the final output includes the nonlinear process detail characteristic data of wave-wave interaction, vortex generation position and energy cascade path, the data format is three-dimensional NetCDF, and the variables include vorticity, strain rate and energy flux. The vertical direction adopts Sigma coordinate layering, a total of 20 layers, and the layer thickness is adaptively adjusted according to the terrain.

[0030] The present application realizes accurate decomposition of nonlinear terms through the cooperation of high-order spectrum method and sub-grid model, separates the nonlinear components of density front and internal wave by using Fourier transform (formula F(k) = ∫u(x)e-ikxdx), solves the problem of insufficient capture of high-frequency oscillation in traditional methods, triggers the sub-grid model through the energy threshold of nonlinear terms (such as 0.001 m 2 / s 2 ), adopts WKB approximation to track phase change, and reveals the energy transfer mechanism from large scale to small scale; combined with the finite volume method to discretize Navier-Stokes equation, the vorticity and energy flux in the front area are calculated, which provides direct observation basis for marine mixing parameterization.

[0031] Step S104, according to the nonlinear process detail characteristic data, the spatio-temporal distribution of density front and internal wave breaking is extracted, the Lagrangian particle tracking method is adopted, the particle is released every 0.1 km, the tracking time step is 1 hour, the mixing coefficient of the front and wave breaking area is calculated, and the fine expression data of marine mixing process is generated.

[0032] Velocity field data is extracted from the current field distribution, and the Lagrangian particle tracking method is used to calculate the particle trajectory through the formula x(t+Δt) = x(t) + u(x, t)Δt, where x is the particle position, u is the velocity field, t is the time, and Δt is the time step, to obtain the particle trajectory distribution. According to the particle trajectory distribution, the velocity gradient field is calculated through the formula The velocity field gradient is calculated, where u is the velocity field and x is the spatial coordinate, to obtain the velocity gradient field distribution.

[0033] Nonlinear features are extracted from the velocity gradient field distribution, and if the absolute value of the velocity gradient exceeds a predetermined threshold, it is marked as a density front area. The characteristic value is calculated through the formula where D is the characteristic value, is the velocity gradient, and the density front area distribution is obtained. According to the density front area distribution, the internal wave breaking feature is calculated, and the spectral analysis method is used to calculate the spectral function through the formula F(k) = ∫u(x)e^(-ikx)dx, where F(k) is the spectral function, u is the velocity field, k is the wave number, and x is the spatial coordinate. If the high-frequency component of the spectral function exceeds a predetermined threshold, it is marked as an internal wave breaking area, and the internal wave breaking area distribution is obtained.

[0034] The spatiotemporal distribution is extracted from the density front area distribution and the internal wave breaking area distribution, and the features are integrated through the formula S(x, t) = D(x, t) + W(x, t), where S is the spatiotemporal distribution, D is the density front feature, W is the internal wave breaking feature, x is the spatial coordinate, and t is the time, to obtain the spatiotemporal distribution data. According to the spatiotemporal distribution data, the regional mixing coefficient is calculated through the formula The mixing coefficient is calculated, where K is the mixing coefficient, is the velocity gradient, and t is the time, to obtain the mixing coefficient distribution.

[0035] Fine expression data is generated from the mixing coefficient distribution, and the grid interpolation method is used to map the mixing coefficient to the grid through the formula G(i, j) = K(xi, yj), where G is the grid feature, K is the mixing coefficient, and xi, yj are the grid coordinates, to obtain the fine expression data of the ocean mixing process.

[0036] In the processing of nonlinear process detail feature data, first, the spatiotemporal distribution data of density front and internal wave breaking are obtained through a high-resolution ocean model, with a spatial resolution of 0.1 kilometers and a time resolution of 1 hour. The Lagrangian particle tracking method is used to release a particle every 0.1 kilometers, with a tracking time step of 1 hour, to ensure that the particle motion trajectory in the density front and wave breaking area is accurately recorded. By calculating the displacement and velocity of the particles, combined with the ocean dynamics equation, the specific location and intensity of the density front and internal wave breaking can be extracted.

[0037] For example, at a certain moment, the displacement of particles in the density front area is 0.5 kilometers, and the speed is 0.1 meters per second, while in the wave breaking area, the displacement is 0.3 kilometers, and the speed is 0.2 meters per second. Based on these data, the mixing coefficient calculation formula is used, such as K = 0.1 * (u^2 + v^2), where u and v are the horizontal and vertical velocities of the particles, respectively, to calculate the mixing coefficients of the front and wave breaking areas.

[0038] For example, in the density front area, the mixing coefficient is 0.01 square meters per second, while in the wave breaking area, the mixing coefficient is 0.02 square meters per second. Finally, these mixing coefficients are combined with the spatiotemporal distribution data to generate fine expression data of the ocean mixing process, providing high-precision basic data support for marine environmental research and prediction.

[0039] Based on Lagrangian particle tracking and velocity gradient analysis, the present application realizes the spatiotemporal visualization of the mixing process, and maps the spatial distribution of density fronts and internal wave breaking (high frequency spectral component ratio > 50%) through particle trajectories (e.g. release interval 0.1 km, tracking time 1 h); based on particle displacement and velocity, the mixing coefficient formula is calculated which solves the limitations of traditional methods relying on empirical formulas; through SVD decomposition, the front and wave breaking characteristics are integrated S(x, t) = D(x, t) + W(x, t)), and a mixing intensity classification map is generated, such as K = 0.01 m 2 / s in the front area and K = 0.02 m 2 / s in the breaking area, realizing the front intensity classification.

[0040] Step S105, according to the fine expression data of the ocean mixing process, the seabed boundary layer turbulence parameterization scheme is optimized, the iterative adjustment method is used to update the turbulent viscosity coefficient in the Mellor-Yamada model, and the optimized turbulent mixing distribution data is generated.

[0041] From the fine expression data, the seabed boundary layer turbulence characteristics are obtained, and the turbulent kinetic energy is calculated by the formula where T is the turbulent kinetic energy, is the velocity gradient, x, y, z are spatial coordinates, and t is time, to obtain the turbulent kinetic energy distribution. According to the turbulent kinetic energy distribution, the initial turbulent viscosity coefficient is calculated using the Mellor-Yamada model, and the formula v = l·sqrt(T) is used, where v is the turbulent viscosity coefficient, l is the mixing length, and T is the turbulent kinetic energy, to obtain the initial turbulent viscosity coefficient distribution.

[0042] If the initial turbulent viscosity coefficient distribution deviates from the turbulent characteristics of the fine expression data by more than a preset threshold, an iterative optimization method is used to update the turbulent viscosity coefficient by the formula v_new = v + a (Tobs - Tsim), where v_new is the updated viscosity coefficient, a is the iteration step, Tobs is the observed turbulent kinetic energy, and Tsim is the simulated turbulent kinetic energy, to obtain the updated turbulent viscosity coefficient distribution.

[0043] According to the updated turbulent viscosity coefficient distribution, the Mellor-Yamada model parameters are adjusted, the turbulent kinetic energy is calculated by the formula q^2 = 2 · T, where q^2 is the turbulent kinetic energy and T is the turbulent kinetic energy, to obtain the adjusted turbulent kinetic energy distribution. The turbulent mixing features are extracted from the adjusted turbulent kinetic energy distribution, and the turbulent mixing intensity is calculated by the formula , where M is the turbulent mixing intensity, v is the turbulent viscosity coefficient, is the velocity gradient, to obtain the turbulent mixing intensity distribution. A grid interpolation method is used to map the turbulent mixing intensity distribution to a regular grid, and the gridded features are calculated by the formula G(i, j, k) = M(xi, yj, zk), where G is the gridded turbulent mixing intensity, M is the turbulent mixing intensity, xi, yj, zk are the grid coordinates, to obtain the optimized turbulent mixing distribution data. If the matching degree of the optimized turbulent mixing distribution data and the fine expression data is lower than a preset threshold, return to step three and repeat the iterative optimization process to obtain the final optimized turbulent mixing distribution data. In the fine expression data based on the ocean mixing process, the least squares optimization algorithm is used to adjust the turbulent viscosity coefficient in the Mellor-Yamada model, and the initial mixing length scale is set to 0.5 meters and the turbulent kinetic energy dissipation rate is set to 1.0 x 10^-6 square meters / second^3. The model parameters are updated by iterative calculation, and the root mean square error of the simulated turbulent field and the observation data is compared in each iteration, and the iteration is stopped when the error is less than 0.001.

[0044] For example, at a depth of 50 meters in the bottom boundary layer, the initial turbulent viscosity coefficient is 0.05 square meters / second, and after 5 iterations, it is optimized to 0.08 square meters / second, and the error is reduced to 0.0008. The K-epsilon turbulent closure scheme is used to calculate the turbulent kinetic energy transport, and the turbulent generation term is set to 1.2 x 10^-5 square meters / second^3 and the dissipation term coefficient is set to 1.92. At a depth of 200 meters

[0045] The turbulent mixing distribution is calculated with a vertical resolution of 0.2 meters within a 200-meter grid area, and the bottom of the boundary layer is processed by a logarithmic law with a roughness length of 0.01 meters. The output of the optimized model shows that in the area of seabed topography mutation, the turbulent mixing coefficient is increased from the background value of 0.01 to 0.03 square meters per second, and the matching degree with the flow velocity profile observed by the ADCP is improved by 15%. The optimized parameters are embedded in the regional ocean model to generate a turbulent mixing distribution dataset covering the entire study sea area with a horizontal resolution of 100 meters and a time interval of 6 hours, which is used to drive high-precision circulation simulation.

[0046] The present application optimizes the turbulent viscosity coefficient from the initial value of 0.05 m 2 / s to 0.08 m 2 / s (error reduction of 15%) through iterative optimization and observation data assimilation based on the Mellor-Yamada model and the RMSE minimization strategy, and realizes dynamic adjustment of the turbulent viscosity coefficient. The introduction of turbulent kinetic energy dissipation rate constraint prevents numerical instability caused by excessive parameterization; through ADCP flow velocity profile comparison, the correlation coefficient of the optimized turbulent mixing coefficient and the measured flow velocity is improved from 0.85 to 0.92, and the root mean square error is reduced to 0.05 m / s, improving the consistency of model-observation.

[0047] Step S106, updating the ocean current field simulation according to the optimized turbulent mixing distribution data, using a three-dimensional non-hydrostatic ocean model to calculate the topography-induced ocean current velocity and vorticity, and generating high-precision ocean current field distribution data.

[0048] The turbulent mixing intensity is extracted from the optimized turbulent mixing distribution data, and calculated by the formula , where M is the turbulent mixing intensity, v is the turbulent viscosity coefficient, is the velocity gradient, x, y, z are spatial coordinates, and t is time, to obtain the turbulent mixing intensity distribution. According to the turbulent mixing intensity distribution, the three-dimensional non-hydrostatic ocean model is initialized, the model equation is discretized by finite volume method, and the initial ocean current velocity is calculated by the formula , where U is the ocean current velocity, p is the pressure, p is the density, and F is the external force term, to obtain the initial ocean current velocity distribution.

[0049] If the matching degree of the initial ocean current velocity distribution and the turbulent mixing intensity distribution is lower than the preset threshold, an iterative adjustment method is used to update the ocean current velocity by the formula Unew = U + β (Mobs - Msim), where Unew is the updated ocean current velocity, β is the adjustment step, Mobs is the observed turbulent mixing intensity, and Msim is the simulated turbulent mixing intensity, to obtain the updated ocean current velocity distribution.

[0050] According to the updated ocean current velocity distribution, the vorticity is calculated by the formula The vorticity distribution is calculated, where ω is the vorticity, U is the ocean current velocity, x, y, z are the spatial coordinates, and t is the time. A grid interpolation method is used to map the ocean current velocity distribution and vorticity distribution onto a regular grid, using the formulas S(i,j,k)=U(xi,yj,zk) and W(i,j,

[0051] The gridded features are calculated using the formula k)=ω(xi,yj,zk), where S is the gridded ocean current velocity, W is the gridded vorticity, and xi, yj, zk are the grid coordinates, thus obtaining the gridded ocean current field distribution.

[0052] If the deviation between the gridded ocean current field distribution and the topographically induced features exceeds a preset threshold, return to step three and repeat the iterative adjustment process to obtain the optimized ocean current field distribution.

[0053] Based on the optimized ocean current field distribution, the topographic induced effect is calculated using the formula. The topographic induced intensity was calculated, where E is the topographic induced intensity, U is the ocean current velocity, h is the topographic height, x, y, z are the spatial coordinates, and t is the time, to obtain high-precision ocean current field distribution data. Specifically, the 3D non-hydrostatic ocean model NEMO was used, with a horizontal grid resolution of 500 meters, 50 vertical σ-coordinate layers, a time step of 30 seconds, and a simulated ocean area covering 100 km × 100 km. The initial ocean current field was driven by optimized turbulent mixed data, using the 3D temperature and salinity field output from the ROMS model as the background field, with an initial flow velocity of 0.1 m / s and an initial vorticity field value of 1.0 × 10^-5 s^-1. The model calculation used a semi-implicit scheme to handle the pressure gradient term, a second-order central difference scheme for momentum equation discretization, and a viscosity coefficient of 0.1 m² / s.

[0054] For regions with varying seafloor topography, a topography-following coordinate transformation was employed. In areas with slopes greater than 5 degrees, the vertical grid was refined to 10 layers, and the horizontal grid was locally refined to 200 meters. During the calculations, the Smagorinsky eddy viscosity parameterization scheme was used to handle sub-grid-scale turbulence, with the eddy viscosity coefficient dynamically adjusted from 0.01 to 0.5 m² / s. After model execution, strong vortex structures were captured in the seafloor canyon region (where water depth drops sharply from 200 meters to 800 meters), with a maximum vorticity of 5.0 × 10⁻⁴ s⁻¹ and a flow velocity increasing to 0.8 m / s. Comparison with HF radar observation data showed a correlation coefficient of 0.92 between the simulation results and measured current directions, with a root mean square error of 0.05 m / s for the velocity. The simulated current field data was stored in NetCDF format at 1-hour intervals, including variables such as three-dimensional velocity, vorticity, and kinetic energy dissipation rate, for subsequent studies on ocean circulation and material transport.

[0055] Step S107, according to the high-precision sea current field distribution data, adjusting the global parameters of the ocean numerical model, including grid resolution and time step, using gradient descent optimization algorithm, if the error of simulation results and observation data is greater than 0.05, adjusting the parameterization scheme, generating improved prediction results data.

[0056] From the high-precision sea current field distribution data, the characteristic parameters are obtained, and the characteristic intensity is calculated by the formula F = Σ | Ui - Umean |, wherein F is the characteristic intensity, Ui is the sea current velocity component, and Umean is the average velocity, and the characteristic parameter distribution is obtained. According to the characteristic parameter distribution, the grid resolution is adjusted by using the grid self-adaptive method, and the updated resolution is calculated by the formula Rnew = Rbase (F / Fmax), wherein Rnew is the updated resolution, Rbase is the basic resolution, F is the characteristic intensity, and Fmax is the maximum characteristic intensity, and the adaptive grid distribution is obtained.

[0057] According to the adaptive grid distribution, the time step is adjusted, and the updated step is calculated by the formula Tnew = Tbase (Rnew / Rmax), wherein Tnew is the updated time step, Tbase is the basic time step, Rnew is the updated resolution, and Rmax is the maximum resolution, and the optimized time step is obtained.

[0058] Using the gradient descent optimization algorithm, according to the optimized time step and the adaptive grid distribution, the simulation error is calculated, and the error value is obtained by the formula J = Σ (Msim - Mobs) ^2, wherein J is the error function, Msim is the simulated sea current field data, and Mobs is the observed sea current field data, and the current simulation error is obtained.

[0059] If the current simulation error is greater than the preset threshold, the parameterization scheme is adjusted, and the model parameters are updated by the formula Pnew = P - α * ΔJ, wherein Pnew is the updated parameter, P is the current parameter, α is the learning rate, is the error gradient, and the updated parameter set is obtained. According to the updated parameter set, the ocean numerical model is run again, and the improved sea current prediction data is generated, and the predicted velocity is calculated by the formula Upred = Usim + ΔU, wherein Upred is the predicted sea current velocity, Usim is the simulated velocity, and ΔU is the velocity correction term, and the improved prediction result data is obtained. According to the improved prediction result data, the matching degree with the observation data is verified, and the verification error is calculated by the formula E = Σ | Upred - Uobs |, wherein E is the verification error, Upred is the predicted sea current velocity, and Uobs is the observed sea current velocity, and the final prediction result data is obtained.

[0060] Based on high-precision sea current field distribution data, the gradient descent optimization algorithm is used to adjust the global parameters of the ocean numerical model. The process is as follows: First, the horizontal grid resolution is gradually optimized from the initial 1000 meters to 300 meters, the vertical layer is adjusted from 30 layers of sigma coordinate to 60 layers, and the time step is shortened from 60 seconds to 20 seconds to improve the spatial and temporal resolution of the model. During the optimization process, the root mean square error between the simulation results and the observation data is calculated. If the error exceeds 0.05, for example, the next grid scale turbulence processing scheme is switched from the Smagorinsky scheme to the Mellor-Yamada scheme, and the dynamic adjustment range of the eddy viscosity coefficient is optimized from 0.02 to 0.4 square meters / second to 0.01 to 0.3 square meters / second. At the same time, for the complex topography of a specific sea area, the adaptive grid refinement technique is used to locally refine the horizontal grid to 100 meters and the vertical grid to 15 layers in the area with a slope greater than 10 degrees. During the optimization process, the simulation error is gradually reduced through iterative calculation, and finally the root mean square error is stabilized within 0.03 meters / second. The output data of the optimized model is stored in NetCDF format with a time interval of 30 minutes, including three-dimensional flow velocity, vorticity, temperature and salinity fields, and turbulent kinetic energy dissipation rate, etc. It provides high-precision data support for the subsequent research on ocean circulation and ecosystem.

[0061] Step S108, according to the improved prediction result data, verify the simulation accuracy of the interaction between complex topography and sea current, use five-fold cross-validation method, compare the fitting degree of simulated sea current field and observation data, generate model precision evaluation data.

[0062] Extract the complex topography features from the prediction result data, calculate the terrain complexity by the formula T = Σ|Zi-Zmean|, where T is the terrain complexity, Zi is the terrain height component, Zmean is the terrain height mean, get the terrain feature distribution. According to the terrain feature distribution, adjust the grid division of the simulated sea current field, calculate the updated grid density by the formula Gnew = Gbase(T / Tmax), where Gnew is the updated grid density, Gbase is the base grid density, T is the terrain complexity, Tmax is the maximum terrain complexity, determine the adaptive grid distribution.

[0063] Using the five-fold cross-validation method, for the simulated sea current field under the adaptive grid distribution, calculate the error of each fold by the formula Ek = Σ|Usim_k-Uobs_k|, where Ek is the error of the kth fold, Usim_k is the simulated sea current velocity of the kth fold, Uobs_k is the observed sea current velocity of the kth fold, get the error data of each fold. Through the random forest algorithm, fuse the error data of each fold and the terrain feature distribution, calculate the weighted features by the formula W = Σw_i F_i, where W is the weighted feature value, w_i is the feature weight, F_i is the terrain feature component, judge the model precision distribution.

[0064] If the variance of the model accuracy distribution is greater than a preset threshold, then the formula is used. Update simulation parameters, where Anew is the updated parameter, A is the current parameter, and β is the adjustment rate. The error gradient is used to obtain the optimized parameter set. Based on the optimized parameter set, simulated ocean current field data is regenerated, and the corrected ocean current velocity is calculated using the formula Upen = Usim + ΔUc, where Upen is the corrected ocean current velocity, Usim is the current simulated velocity, and ΔUc is a correction term based on terrain features, generating improved ocean current field data. For the improved ocean current field data, a five-fold cross-validation method is used, using the formula Efinal = ...

[0065] The final validation error is calculated using Σ|Unew-Uobs|, where Efinal is the final validation error, Unew is the corrected ocean current velocity, and Uobs is the observed ocean current velocity, thus obtaining the final model accuracy evaluation data.

[0066] To verify the accuracy of the simulation of the interaction between complex terrain and ocean currents, a five-fold cross-validation method was employed, comparing and analyzing the improved prediction results with observational data. First, the simulated ocean current field data and observational data were divided into five subsets according to time series, each containing data for 30 consecutive days to ensure temporal continuity. For each subset, the root mean square error (RMSE) of the simulation results and observational data was calculated, and the error distribution was statistically analyzed.

[0067] For example, in the first subset, the RMSE of the simulated ocean current field was 0.04 m / s, while in the second subset, the RMSE decreased to 0.035 m / s, indicating that the simulation accuracy of the model varied across different time periods. To further analyze the sources of error, Principal Component Analysis (PCA) was used to reduce the dimensionality of the simulation results and observation data, extracting the main feature vectors. It was found that the error contribution rate was high in complex terrain areas, reaching 45%. To address this phenomenon, a terrain correction factor was introduced, and the local grid was refined in areas with a terrain slope greater than 15 degrees, increasing the horizontal resolution from 200 meters to 50 meters and the vertical resolution from 20 layers to 30 layers. Simultaneously, the Kriging interpolation algorithm was used to spatially interpolate the observation data, generating a high-resolution observation field, which was then compared point-by-point with the simulation results. Through iterative calculations, the model parameters were gradually optimized, ultimately stabilizing the average RMSE of the five-fold cross-validation to within 0.03 m / s, significantly improving the model accuracy. The optimized model accuracy evaluation data is stored in CSV format, containing fields such as time series, spatial location, simulated values, observed values, and error values, providing reliable data support for subsequent research.

[0068] Step S109, according to the model precision evaluation data, feedback optimization adaptive grid partitioning technology and turbulence closure model, if the precision is lower than 0.9, adjust the grid density to the minimum unit 0.05 kilometers and Mellor-Yamada model turbulence viscosity coefficient, generate the final optimization simulation scheme data.

[0069] By analyzing the model precision evaluation data, the precision distribution characteristics are calculated, and the precision deviation is calculated by the formula P = Σ | Di-Dmean |, wherein P is the precision deviation, Di is the single-point precision value, and Dmean is the precision mean value, to obtain the precision distribution characteristics. According to the precision distribution characteristics, it is judged whether the model precision is lower than the preset threshold value, if lower than the preset threshold value, the grid density is adjusted by the formula Gadj = Ginit (P / Pmax), wherein Gadj is the adjusted grid density, Ginit is the initial grid density, P is the precision deviation, and Pmax is the maximum precision deviation, to determine the adjusted grid density. According to the adjusted grid density, the adaptive grid partitioning is updated, the grid partitioning parameter is calculated by the formula Mnew = Mbase + a (Gadj-Ginit), wherein Mnew is the updated grid parameter, Mbase is the basic grid parameter, a is the adjustment coefficient, Gadj is the adjusted grid density, and Ginit is the initial grid density, to obtain the updated grid partitioning parameter. The Mellor-Yamada turbulence model is optimized by the updated grid partitioning parameter, and the turbulence viscosity coefficient is adjusted by the formula Vnew = Vinit (1+βP), wherein Vnew is the adjusted viscosity coefficient, Vinit is the initial viscosity coefficient, β is the adjustment factor, and P is the precision deviation, to determine the optimized turbulence model parameter.

[0070] According to the optimized turbulence model parameter, the simulation sea current field data is generated, the optimized sea current velocity is calculated by the formula Uopt = Ubase + γ (Vnew-Vinit), wherein Uopt is the optimized sea current velocity, Ubase is the basic sea current velocity, γ is the correction factor, Vnew is the adjusted viscosity coefficient, and Vinit is the initial viscosity coefficient, to obtain the optimized sea current field data. According to the optimized sea current field data, the K-means algorithm clustering analysis is adopted to generate the precision distribution cluster, and the cluster error is calculated by the formula C = Σ | Uopt_k-Uobs_k |, wherein C is the cluster error, Uopt_k is the kth cluster optimized sea current velocity, and Uobs_k is the kth cluster observed sea current velocity, to determine the precision distribution cluster characteristics. According to the precision distribution cluster characteristics, the final optimization simulation scheme data is generated, the scheme parameter is adjusted by the formula Sfinal = Sinit + δC, wherein Sfinal is the final scheme parameter, Sinit is the initial scheme parameter, δ is the update factor, and C is the cluster error, to obtain the final optimization simulation scheme data.

[0071] Based on the model accuracy evaluation data, when the accuracy is less than 0.9, the system automatically triggers the adaptive grid optimization process, adopts the Delaunay triangulation algorithm to dynamically encrypt the calculation domain, and gradually refines the grid size from the default 0.1 km to 0.05 km in the area where the flow velocity gradient is greater than 0.1 s-1, while adopting σ coordinate transformation along the vertical direction to increase the number of layers to 40 layers. For the optimization of the turbulence closure model, the Mellor-Yamada model parameters are inverted by the Levenberg-Marquardt algorithm, and the initial value of the turbulent viscosity coefficient is set to 0.2 m 2 / s, and the deviation between the surface 10 m flow velocity observation data and the simulation value is dynamically adjusted, when the surface RMSE exceeds 0.05 m / s, the viscosity coefficient is increased by 0.01 m 2 / s step, and the upper limit is set to 0.5 m 2 / s. The shear production item correction coefficient is introduced in the eddy energy equation, and the velocity gradient tensor is calculated by the second-order central difference, when the strain rate is greater than 0.05 s-1, the coefficient is adjusted from 1.0 to 1.2. During the optimization process, the Richardson number distribution is monitored in real time, and double time step calculation is implemented in the mixed layer area where Ri<0.25. The final output includes the encrypted grid topology, the corrected turbulent parameter field and the three-dimensional NetCDF data set of the space-time distribution of each physical quantity, and the spatial correlation coefficient of the horizontal flow velocity field is improved to 0.92, and the average error of the vertical temperature profile is controlled within 0.15℃.

[0072] The above-described embodiments are merely used to illustrate the technical solutions of the present application, but not limit the same; even though the present application has been described in detail with reference to the foregoing embodiments, those ordinarily skilled in the art should understand: the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for simulating the interaction between seabed topography and ocean currents, characterized in that, include: Acquire high-resolution seabed topography data, generate an adaptive grid structure, and construct refined topography description data; Based on the refined terrain description data, a turbulent closure model is used to generate interaction data between the ocean current field and the terrain; Based on the interaction data, a nonlinear process is simulated to generate nonlinear process characteristic data; Based on the nonlinear process characteristic data, ocean mixing process data are generated using a particle tracking method. Based on the ocean mixing process data, turbulence parameters are optimized to generate optimized turbulence mixing data; Based on the optimized turbulent mixing data, the ocean current field simulation is updated to generate high-precision ocean current field distribution data; The acquisition of high-resolution seabed topographic data, generation of an adaptive grid structure, and construction of refined topographic description data include: Obtain a set of three-dimensional coordinate points from a high-resolution seabed topography dataset; Calculate the terrain gradient based on the set of three-dimensional coordinate points; If the terrain gradient is greater than a preset threshold, the quadtree algorithm is used to encrypt the mesh. Through the encrypted grid, a refined terrain description data is generated using the bilinear interpolation method; The step of simulating a nonlinear process and generating nonlinear process characteristic data based on the interaction data includes: Obtain a nonlinear feature distribution from the interaction data; The energy of the nonlinear term is calculated using a higher-order spectral method. If the energy of the nonlinear term is greater than a preset threshold, a subgrid model is introduced; Using the sub-mesh model and the finite volume method, nonlinear process characteristic data are generated. The step of optimizing turbulence parameters and generating optimized turbulence mixing data based on the ocean mixing process data includes: Turbulence characteristics are obtained from the ocean mixing process data; The initial turbulent viscosity coefficient is calculated using a closed turbulence model. If the initial turbulent viscosity coefficient deviates from the turbulent characteristics by a greater than a preset threshold, an iterative optimization method is used to update the turbulent viscosity coefficient and generate optimized turbulent mixing data.

2. The method as described in claim 1, characterized in that, The step of generating interaction data between ocean currents and topography based on the refined topographic description data and using a turbulent closure model includes: Extract a set of terrain feature points from the refined terrain description data; The turbulence mixing coefficient is calculated using a closed turbulence model. Extract ocean current field shear characteristics from the turbulent mixing coefficient; Principal component analysis was used to reduce the dimensionality of the interactive feature vectors and generate data on the interaction between ocean currents and topography.

3. The method as described in claim 1, characterized in that, The step of generating ocean mixing process data using a particle tracking method based on the nonlinear process characteristic data includes: Extract velocity field data from the nonlinear process characteristic data; The Lagrange particle tracking method is used to calculate the particle trajectory distribution; Calculate the velocity gradient field based on the particle trajectory distribution; Density fronts and internal wave breaking features are extracted from the velocity gradient field to generate ocean mixing process data.

4. The method as described in claim 1, characterized in that, The iterative optimization method is used to update the turbulent viscosity coefficient and generate optimized turbulent mixing data, including: Gradient descent optimization strategy: The gradient descent algorithm is adopted, with the turbulent viscosity coefficient as the optimization variable and the root mean square error between the simulated turbulent field and the observed data as the loss function to calculate the gradient direction; Adjust terrain adaptive step size: dynamically adjust and optimize step size according to terrain gradient. If the terrain gradient is greater than the preset threshold, the step size will be reduced to 50% of the original step size to prevent parameter oscillation in steep terrain areas. Controlling Iteration Convergence: Set the maximum number of iterations and the convergence threshold, repeat the iteration until the conditions are met, and generate the final optimized turbulent mixing data.

5. The method as described in claim 4, characterized in that, The iterative optimization method further includes: Dynamic coupling of multiple parameters: Based on gradient descent optimization, turbulent kinetic energy dissipation rate is introduced as an auxiliary optimization variable to establish the nonlinear relationship between turbulent viscosity coefficient and . Parameter sensitivity screening: Five-fold cross-validation is used to screen parameter combinations that are sensitive to simulation errors, generating robust turbulent mixing data.

6. A simulation system for the interaction between seabed topography and ocean currents, implementing the simulation method of claim 1, comprising: The terrain acquisition unit is used to acquire high-resolution seabed terrain data, generate an adaptive grid structure, and construct refined terrain description data. An interaction data generation unit is used to generate interaction data between ocean currents and topography based on the refined topographic description data and a turbulent closure model. The nonlinear process feature data generation unit simulates the nonlinear process based on the interaction data and generates nonlinear process feature data. The ocean mixing process data generation unit generates ocean mixing process data using a particle tracking method based on the nonlinear process characteristic data. An optimized turbulent mixing data generation unit optimizes turbulence parameters and generates optimized turbulent mixing data based on the ocean mixing process data. The ocean current field distribution data generation unit updates the ocean current field simulation based on the optimized turbulent mixing data, and generates high-precision ocean current field distribution data.

7. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed, implements the simulation method for the interaction between seabed topography and ocean currents as described in any one of claims 1-5.

Images