A numerical method for simulating tsunami and lake surge based on Boussinesq-type equation

Abstract

The application discloses a numerical simulation method for tsunamis and lake surges based on Boussinesq type equations, and comprises the following steps: obtaining digital elevation and initial disturbance data of a target area; constructing a high-order Boussinesq type control equation containing adjustable dispersion parameters and a frequency dispersion term to adjust the dispersion approximation precision; under a unified explicit time advancing system, discretizing the conservation term and the frequency dispersion term by using a finite volume method and a finite difference method respectively, and loading data to a calculation grid; judging the dry and wet states based on a preset water depth threshold, and truncating the frequency dispersion term for solving in the dry grid area to obtain simulation results of a flow field and a water level field. By constructing the control equation containing the adjustable dispersion parameters and the mixed numerical discrete framework, combining the explicit time advancing and the dry and wet grid truncation processing, the method does not need iteration and is simple to calculate; the dispersion description precision is improved, and the numerical stability problem of the dry and wet boundary is solved.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology of water wave dynamics, and in particular to a numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations. Background Technology

[0002] Nearshore water waves exhibit strong nonlinear and dispersion characteristics. Traditional shallow water equations (SWEs) based on the hydrostatic pressure assumption fail to reflect the dispersion effects during wave propagation due to neglecting the influence of vertical acceleration. They also struggle to accurately describe the differences in propagation velocity of different wavelength components, thus limiting their application in simulating wave-like tidal bores and secondary wave train structures. To overcome the hydrostatic pressure limitation, related research has introduced Boussinesq terms to describe dispersion effects by considering higher-order influences on the vertical velocity distribution of water waves, forming Boussinesq-type equations. Among them, Beji and Nadaoka proposed a two-dimensional higher-order Boussinesq-type equation applicable to complex terrain conditions based on the classical Boussinesq equations. By modifying the nonlinear and dispersion terms, the equations can describe the effects of non-hydrostatic pressure and wave dispersion effects, and expand the applicable water depth range of the model, featuring accurate dispersion calculation and a wider applicability in deep water. In terms of numerical solutions, Soares-Frazão and Zech proposed a hybrid numerical discretization method, which uses the finite volume method to discretize the conservation terms in the governing equations to ensure conservation, and uses the finite difference method to handle higher-order dispersion terms to flexibly deal with higher-order derivatives, thereby realizing the joint solution of nonlinear and dispersion terms.

[0003] Numerical models of water waves based on Boussinesq theory have been implemented in various forms internationally. For example, the FUNWAVE model, based on fully nonlinear Boussinesq equations, is used to simulate nearshore wave propagation and strong dispersion processes. The JAGURS-D model introduces Boussinesq dispersion corrections on the basis of traditional models for high-precision tsunami numerical simulation. In earlier technical solutions, water wave and tsunami simulations mainly adopted models based on shallow water wave equations, such as TUNAMI, COMCOT, and Clawpack. These methods are based on hydrostatic pressure assumptions and are mainly used for large-scale long-wave propagation simulations, but they have certain limitations in describing complex nearshore water wave structures.

[0004] However, existing technologies still have the following shortcomings in numerical simulation of water waves: In terms of computational complexity, although some existing Boussinesq-type models can consider nonlinear and dispersion effects, the higher-order derivative terms in their governing equations often require the introduction of complex numerical discretization schemes or multi-level iterative solution processes, resulting in computational difficulties and low computational efficiency, which is not conducive to engineering applications of large-scale or long-term processes; In terms of describing dispersion characteristics, some existing Boussinesq models still have shortcomings in the accuracy of approximating dispersion relations, and there is a deviation between their linear dispersion relations and theoretical solutions, affecting the accurate characterization of actual water wave evolution processes; In terms of handling dry and wet boundaries, existing technologies often have difficulty effectively handling the dynamic changes of dry and wet grids when simulating the propagation and receding of water bodies involved in tidal bore propagation. Some models are prone to numerical instability, non-physical oscillations, or mass non-conservation phenomena in the dry-wet boundary region, which is particularly prominent under complex estuarine topography or strongly nonlinear flow conditions. Therefore, there is a need for a numerical simulation method that can improve the accuracy of describing dispersion effects while ensuring computational simplicity, and achieve stable handling of dynamic changes in dry and wet grids, so as to perform high-precision and high-stability numerical simulations of complex nonlinear water wave processes such as tidal bores, dam breaks, and tsunamis. Summary of the Invention

[0005] The purpose of this invention is to provide a numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations. By constructing control equations with adjustable dispersion parameters and establishing a hybrid numerical discretization framework, combined with a unified explicit time-progression system and dynamic truncation processing of wet and dry grids, the method ensures computational simplicity without the need for iterative solutions, while improving the dispersion description accuracy of propagation characteristics of different wavelength components and solving the numerical stability problem in dynamically changing regions of wet and dry boundaries.

[0006] To address the aforementioned technical problems, a first aspect of this invention provides a numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations, comprising the following steps: Step S1: Obtain digital elevation data and initial disturbance source data for the target area; Step S2: Construct a control equation containing an adjustable dispersion parameter. The control equation is a high-order Boussinesq type equation that introduces a dispersion term. The adjustable dispersion parameter is used to adjust the approximation accuracy of the control equation for the difference in propagation speed of different wavelength components. Step S3: Establish a hybrid numerical discretization framework, load the digital elevation data and the initial disturbance source data into the computational grid corresponding to the control equation, and discretize the conservation terms in the control equation using the finite volume method and the dispersion terms in the control equation using the finite difference method under a unified explicit time-progression system. Step S4: During the calculation process of the hybrid numerical discretization framework, a dry-wet grid dynamic processing is performed to obtain the numerical simulation results of the flow field and water level field. The dry-wet grid dynamic processing includes: determining the dry-wet state of the grid cell based on a preset water depth threshold, and truncating the solution of the dispersion term in the region involving the dry grid cell.

[0007] Furthermore, step S3 involves discretizing the dispersion term using the finite difference method, including: Step S331: Split the spatiotemporal mixed partial derivatives contained in the dispersion term, and transform the mixed term of the second-order time derivative and spatial derivative into an algebraic combination of the first-order spatial derivative and the second-order spatial derivative. Step S332: On the computational grid, differentially discretize the higher-order spatial derivatives in the dispersion term and construct a local semi-implicit dispersion correction expression; Step S333: Within the current time step, substitute the computable expression as a known source term into the governing equation to calculate the value of the dispersion term.

[0008] Furthermore, the adjustable dispersion parameter participates in the construction of the computable expression in step S333 in the form of a scaling factor; When the adjustable dispersion parameter takes different values, the contribution strength of the dispersion term in the control equation is adjusted accordingly, so that the control equation switches between the Peregrine dispersion form and the improved dispersion relationship model.

[0009] Further, in step S3, loading the initial disturbance source data into the computational grid corresponding to the control equation includes: Step S311: The simulation framework is composed of the governing equations, the hybrid numerical discretization framework, and the dry and wet mesh dynamic processing. The solver structure of the simulation framework remains unchanged in the following steps. Step S312: Select the corresponding initial disturbance source data according to the type of the target water wave scene; wherein, the earthquake tsunami scene corresponds to the earthquake fault parameter data, the tidal bore scene corresponds to the time-varying boundary data generated by the boundary condition function, and the landslide lake tsunami scene corresponds to the landslide motion parameter data. Step S313: Load the selected initial disturbance source data into the computational grid; Step S314: When the target water wave scene is switched, replace the type of the initial disturbance source data and the corresponding source item processing operation, and return to execute step S312.

[0010] Furthermore, the initial disturbance source data is seismic fault parameter data, which is stored in a data file format independent of the solver code; Step S313, loading the initial disturbance source data into the computational grid corresponding to the control equation, includes: Step S3131: Read the data file, obtain the wave source position parameters specified in the data file in the form of grid coordinates, and map the wave source position to the computation grid; Step S3132: When the resolution of the computational grid changes, the geographical location of the wave source in the target area is kept unchanged by adjusting the grid coordinate parameters in the data file.

[0011] Furthermore, the initial disturbance source data is earthquake fault parameter data, which includes: source location, fault strike, dip angle, slip angle, fault length, fault width, and / or slip amount.

[0012] Furthermore, in step S3, the conserved terms are discretized using the finite volume method, including: Step S321: Extract the physical quantities of the mesh cells on both sides of the interface, and calculate the local gradient ratio of the physical quantities; Step S322: Input the local gradient ratio into the flux limiter function, and the flux limiter function determines the interface interpolation weights used to calculate the flux; Step S323: Based on the interface interpolation weight, the physical quantities on both sides of the interface are reconstructed by weighting to obtain the interface reconstruction value used to calculate the flux; Step S324: Substitute the interface reconstruction value into the flux calculation formula to calculate the flux through each grid interface, and use the flux to update the conserved quantity.

[0013] Furthermore, the flux limiter function is in upwind, QUICK, or Superbee format; The flux limiter function uses higher-order interpolation weights in low-gradient regions and lower-order interpolation weights in high-gradient regions, based on the magnitude of the local gradient ratio.

[0014] Further, obtaining the digital elevation data in step S1 includes: Step S111: Obtain digital elevation data of the study area in its original format; Step S112: Convert the original format digital elevation data into a three-column data file of a unified format, wherein the three-column data file sequentially includes a longitude data column, a latitude data column, and a topographic elevation data column; Step S113: Using the three-column data file as the input format for the digital elevation data, load it into the computational grid.

[0015] Furthermore, the step S4, which involves truncating the solution of the dispersion term in the region involving dry mesh cells, includes: Step S41: Before performing the differential solution of the dispersion term on the current wet mesh cell, determine all adjacent mesh cells covered by the dispersion term differential template; Step S42: Detect the dry / wet status indicators of all adjacent grid cells one by one; Step S43: If any dry grid cell exists among all adjacent grid cells, set the calculated value of the dispersion term at the current wet grid cell to zero and skip the subsequent solution steps for the dispersion term. Step S44: If all adjacent mesh cells are wet mesh cells, then calculate according to the differential combination expression of the dispersion term.

[0016] Accordingly, a third aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described numerical simulation method for tsunamis and lake tsunamis based on Boussinesq type equations.

[0017] Accordingly, a fourth aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations.

[0018] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By splitting the spatiotemporal mixed partial derivatives in the dispersion term and transforming them into an algebraic combination of the first and second spatial derivatives, this invention constructs a locally semi-implicit time-propelling architecture. This allows higher-order dispersion terms to be directly substituted into the governing equations as known source terms. At the same time, adjustable dispersion parameters are introduced to participate in the dispersion term calculation as scaling factors. Thus, under a unified explicit time-propelling framework, the approximation accuracy of the dispersion relationship can be dynamically adjusted without iterative solutions. This achieves an accurate description of the propagation characteristics of different wavelength components while significantly reducing computational complexity. 2. By introducing a dynamic identification mechanism for wet and dry grids based on a preset water depth threshold into the hybrid numerical discretization framework, and by checking the state of adjacent grid cells covered by the differential template one by one before calculating the dispersion term, the solution of the dispersion term is truncated once a dry grid cell is found. This effectively curbs the non-physical numerical oscillations caused by data loss in the wet and dry interface area, and ensures the stability and reliability of the simulation results of the flow field and water level field in the whole computational domain during the advance and fall of the water body. 3. By constructing a simulation framework composed of governing equations, a hybrid numerical discretization framework, and dynamic processing of wet and dry grids, and storing earthquake fault parameters in independent data files or defining initial disturbance sources for different scenarios using boundary condition functions and landslide motion parameters, the same solver structure can maintain the geographical location of the wave source at different grid resolutions and flexibly switch between different water wave scenarios. Combined with the flux limiter reconstruction of interface flux in the finite volume method discretization, numerical overshoot and dissipation in high gradient regions are further suppressed. Overall, high-precision and robust engineering simulation of various complex hydrodynamic processes such as tsunamis, tidal bores, and landslide tsunamis is achieved. Attached Figure Description

[0019] Figure 1 This is a flowchart of the numerical simulation method for tsunamis and lake tsunamis based on Boussinesq type equations provided in this embodiment of the invention; Figure 2 This is a schematic diagram of the free surface evolution (t = 0 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 3 This is a schematic diagram of the free surface evolution (t = 1 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 4 This is a schematic diagram of the free surface evolution (t = 20 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 5 This is a schematic diagram of the free surface evolution (t = 40 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 6 This is a schematic diagram of the free surface evolution (t = 60 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 7 This is a schematic diagram of the free surface evolution (t = 80 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 8 This is a schematic diagram of the free surface evolution (t = 100 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 9 This is a schematic diagram of the free surface evolution (t = 120 min) during the propagation of the Qiantang River tidal bore, provided in an embodiment of the present invention. Figure 10 This is a schematic diagram of the free surface evolution (t = 0 s) of landslide-induced waves provided in an embodiment of the present invention; Figure 11 This is a schematic diagram of the free surface evolution (t = 13 s) of landslide-induced waves provided in an embodiment of the present invention; Figure 12 This is a schematic diagram of the free surface evolution (t = 25 s) of landslide-induced waves provided in an embodiment of the present invention; Figure 13 This is a schematic diagram of the free surface evolution (t = 35 s) of landslide-induced waves provided in an embodiment of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0021] Please refer to Figure 1 The first aspect of this invention provides a numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations, comprising the following steps: Step S1: Obtain digital elevation data and initial disturbance source data for the target area.

[0022] When implementing the water wave numerical simulation method provided in this embodiment, it is first necessary to prepare basic input data for subsequent mathematical modeling and numerical solution, namely, to obtain digital elevation data of the target study area and initial disturbance source data corresponding to the type of water wave problem in that area. Digital elevation data describes the spatial distribution of seabed topography or riverbed elevation in the target area, and is usually derived from global or regional topographic databases, such as global topographic relief data with a resolution of 15 arcminutes. To unify the data interface, the raw format digital elevation data can be converted into a standard three-column text file arranged in latitude coil order, with each line containing longitude, latitude, and topographic elevation values, enabling the program to accurately read the data according to the defined grid structure. Initial disturbance source data is used to define the original dynamics that trigger water wave motion. Its physical meaning depends on the specific scenario to be simulated: when simulating a tsunami triggered by an undersea earthquake, the data may be a parameter file describing the fault rupture characteristics, which at least implicitly contains information such as the location of the epicenter, fault geometry, and slip amount. These parameters will be used to construct the initial water surface disturbance within the computational domain; when simulating tidal bore propagation, the data can reflect the time-varying function controlling the boundary water level changes, used to drive water movement at the boundary of the computational domain; when simulating a lake tsunami caused by a landslide entering the water, the data can reflect the mass distribution and kinematic behavior of the landslide. By standardizing both digital elevation data and initial disturbance source data during the preparation stage, the consistency and repeatability of the computational configuration in subsequent steps are ensured.

[0023] Step S2: Construct a control equation containing an adjustable dispersion parameter. The control equation is a high-order Boussinesq type equation that introduces a dispersion term. The adjustable dispersion parameter is used to adjust the accuracy of the control equation in approximating the differences in propagation speeds of different wavelength components.

[0024] After acquiring and standardizing the input data, it is necessary to establish governing equations describing the physical laws of water wave motion. The constructed governing equations are a class of high-order Boussinesq-type equations, which introduce a dispersion term to characterize the effects of non-hydrostatic pressure on the basis of traditional shallow water equations, overcoming the limitation that the hydrostatic pressure assumption cannot reflect the influence of vertical acceleration. Mathematically, this dispersion term is a mixed high-order partial derivative term containing the second-order time derivative and spatial derivative. Its existence enables the equations to describe the physical fact that water waves of different wavelengths propagate at different speeds, i.e., the dispersion effect. To reduce the approximation deviation between the model's dispersion relation and the exact linear wave theory, the governing equations include an adjustable dispersion parameter. This parameter is not derived from a specific physical constant but is embedded in the construction of the dispersion term as a dimensionless coefficient. By adjusting the value of this parameter, the dispersion relation curve after linearization of the equations can be changed, allowing it to smoothly transition from the classic Peregrine-type dispersion form to an improved model that more precisely approximates the complete Airy wave dispersion relation. This adjustable mechanism avoids the need to re-derive equations for different water depths or wavelength ranges, giving a single governing equation framework the ability to maintain high descriptive accuracy under various water wave dispersion characteristics. In practical applications, the core equation structure can be kept unchanged while only the value of the parameter is changed, allowing the same solution program to adapt to different simulation needs, from deep-sea tsunami propagation to nearshore strong dispersion tidal evolution.

[0025] Step S3: Establish a hybrid numerical discretization framework, load digital elevation data and initial disturbance source data into the computational grid corresponding to the control equations, and discretize the conservation terms in the control equations using the finite volume method and the dispersion terms using the finite difference method under a unified explicit time-progression system.

[0026] First, the digital elevation data obtained in step S1 is interpolated onto a predefined structured computational grid to form a static topographic field. Simultaneously, the initial disturbance source data is mapped to the corresponding grid locations according to the simulated physical process, generating the initial free surface disturbance or velocity field distribution, thus providing initial conditions for the time-progression calculation. For spatial discretization, the framework employs a hybrid strategy combining the finite volume method and the finite difference method: for convection and flux terms representing mass and momentum conservation in the governing equations, the finite volume method is used to calculate numerical fluxes at the grid cell interfaces, ensuring strict conservation of water volume and momentum during the calculation process; for dispersion terms containing complex derivative structures, the finite difference method is used to directly perform spatial partial derivative discretization approximations at adjacent grid nodes, thus flexibly handling higher-order terms that are difficult to express uniformly using the finite volume method. In terms of time progression, the entire solution process is built within a unified explicit framework. Within each time step, the intermediate water level field is first calculated based on the continuity equation and the velocity field at the current moment. Then, the updated water level field is used to sequentially evaluate the physical contributions of the pressure gradient term, convection term, and dispersion term in the momentum equation, and the velocity field at the new moment is directly obtained by explicit superposition. Finally, the continuity equation is used again to complete the final update of the water level. This multi-step update mechanism completely avoids the process of constructing the overall stiffness matrix or iteratively solving linear equations, making the program implementation simple and easy to extend on parallel computing architectures.

[0027] Step S4: During the calculation process of the hybrid numerical discretization framework, perform dynamic processing of dry and wet grids to obtain the numerical simulation results of the flow field and water level field. The dynamic processing of dry and wet grids includes: determining the dry and wet state of the grid cells based on a preset water depth threshold, and truncating the solution of the dispersion term in the region involving dry grid cells.

[0028] During the advance and fall of water, the area traversed by the wavefront exhibits a phenomenon where grid cells alternate above and below the water surface, forming a dynamic wet-dry boundary. If hydrodynamic terms are calculated directly in the wet-dry transition zone without processing, especially for the difference operation of dispersion terms involving multiple adjacent grid values, non-physical numerical oscillations can easily occur due to missing data. To address this issue, a dynamic wet-dry grid processing mechanism is embedded in the computation process of the hybrid numerical discretization framework. This mechanism sets a water depth threshold for each grid cell and, after updating the water level at each time step, checks in real-time whether the water depth of each grid cell is greater than the threshold. If the water depth is lower than the threshold, the cell is marked as a dry grid, and its water depth is forcibly set to zero, preventing it from participating in subsequent momentum exchange calculations. More importantly, when performing the difference solution of dispersion terms on a wet grid, the state of all adjacent grid cells covered by its difference template must be checked beforehand. Once any dry grid cell is detected in an adjacent cell, the dispersion term calculation at that wet grid cell is truncated. Specifically, the contribution of the dispersion term at that location in that time step is directly set to zero, and the corresponding difference summation step is skipped. The normal dispersion term solution process only proceeds when all relevant adjacent cells are wet grid cells. This approach of setting a boundary protection threshold at the most computationally intensive high-order dispersion term stage effectively prevents spurious oscillations introduced by incomplete dry / wet boundary data from propagating into the computational domain. This ensures that the simulation of the entire flow field and water level field can proceed smoothly under complex shoreline and dynamic boundary conditions, ultimately yielding physically reasonable numerical simulation results for the flow field and water level field.

[0029] Through the implementation of the above steps, this invention, within a unified explicit time-progressive hybrid discrete framework, allows for flexible adjustment of the dispersion relation approximation accuracy of the control equations using adjustable dispersion parameters. Simultaneously, it employs the finite volume method and the finite difference method to handle the conservation and dispersion terms respectively, and conditionally truncates the dispersion term in the wet-dry boundary region. This simplifies the entire solution process without requiring iteration or additional storage variables, significantly improving the accuracy of describing the wavelength-dependent propagation characteristics of complex water waves. Furthermore, it effectively solves the problem of instability in the calculation of wet-dry boundaries during water level fluctuations. Therefore, it can provide high-precision and stable flow field and time-varying water level simulation results for various unsteady water wave processes such as tidal bores, tsunamis, and landslide tsunamis.

[0030] Specifically, step S1 involves acquiring digital elevation data, including: Step S111: Obtain digital elevation data of the study area in its original format.

[0031] In practical implementation, the process of acquiring digital elevation data in step S1 can be further refined into a complete data preprocessing and standardized loading process. First, the original format of the digital elevation data for the study area is acquired. The source and format of this original data vary depending on the application scenario: In global tsunami simulations, gridded topographic relief files, such as global topographic data with a resolution of 15 or 10 arcminutes, can be downloaded from a global topographic database. The original storage format is typically the gridded binary format commonly used for geospatial data. In estuary tsunami simulations, the original data may originate from measured scatter data or local high-resolution digital elevation models provided by waterway surveying departments. In reservoir landslide tsunami simulations, underwater topographic survey results or hydraulic engineering survey data may be used. These original data differ in coordinate system, data arrangement, and file format. If they are directly read and parsed by the solver, it not only increases program complexity but also easily introduces data mapping errors due to format incompatibility.

[0032] Step S112: Convert the original format digital elevation data into a unified format three-column data file, which contains longitude data column, latitude data column and terrain elevation data column in sequence.

[0033] After acquiring the raw data, a format standardization process is performed to convert digital elevation data from different sources into a unified three-column data file. This three-column file is a plain text structure, with each row corresponding to a terrain sampling point. It contains columns for longitude, latitude, and terrain elevation, separated by spaces or tabs. The entire set of sampling points is arranged row by row according to latitude circles from north to south. This latitude-longitude priority arrangement aligns with the subsequent practice of indexing the computational grid by longitude and latitude, facilitating rapid retrieval of terrain elevation values ​​at any node using the grid index. During the format conversion process, if the resolution of the raw data is higher than the target computational grid resolution, local averaging or bilinear interpolation can be used to merge the elevation information of densely sampled points to the computational node locations, avoiding wasting data precision. If the raw data contains gaps or outliers in certain areas, gap filling or outlier removal can be performed during conversion to ensure that the final terrain field loaded into the computational grid has no abrupt breaks. For areas bounded by coastlines, land and water can be distinguished by specific elevation marker values ​​in three data files, providing initial topographical basis for subsequent dry and wet grid determination.

[0034] Step S113: Load the three-column data file as the input format for digital elevation data into the computational grid.

[0035] After format conversion, the three-column data file is loaded into the computational grid as a standardized input format for digital elevation data. During loading, the program reads the file once, statistically analyzing the sampling range and step size in the longitude and latitude directions, and accordingly creating a two-dimensional elevation array in memory that corresponds one-to-one with the computational grid. During the second read, the program searches for the elevation value of the nearest sampling point from the three-column file based on the longitude and latitude coordinates of each grid node and fills the corresponding position in the array. For grid nodes located at the edge of the computational domain or outside the sampling points, reasonable elevation values ​​can be assigned through extrapolation or constant filling. At this point, the digital elevation data is loaded in the form of a structured array, which can be directly used in step S3 to construct the topographic field and for subsequent physical calculations. This unified format preprocessing and standardized loading process allows raw topographic data from different sources and resolutions to be accessed through the same interface into the hybrid numerical discretization framework, ensuring not only the accuracy and consistency of topographic input but also providing standardized data support for flexible switching between different water wave scenarios.

[0036] Specifically, step S3 involves loading the initial disturbance source data into the computational grid corresponding to the governing equations, including: Step S311: The simulation framework consists of the governing equations, the hybrid numerical discretization framework, and the dynamic processing of wet and dry meshes. The main structure of the solver of the simulation framework remains unchanged in the following steps.

[0037] Before loading the initial disturbance source data, the main simulation framework was constructed and solidified. This framework consists of three core functional modules: the governing equation module defines a mathematical model based on high-order Boussinesq equations, responsible for describing the nonlinear and dispersion effects in the water wave propagation process; the hybrid numerical discretization framework module provides spatial and temporal discretization strategies, including finite volume method discretization of conservation terms and finite difference method discretization of dispersion terms, as well as a unified explicit time-progression mechanism; and the wet-dry mesh dynamic processing module undertakes the numerical stability control task of the moving boundary region. These three parts work together to form a complete numerical solution system, establishing a full-link technical solution from physical equations to numerical algorithms to boundary processing. Unlike conventional practices that develop dedicated solvers for different water wave problems, this framework solidifies the main solver structure into a general computational kernel independent of the specific disturbance cause. This main structure remains stable during subsequent initial disturbance loading and scene switching operations for different scenarios. This design allows the core algorithm components in the solver's main structure, such as the grid data structure, discretization format selection, time step progression logic, and wet / dry determination rules, to be adapted to different water wave simulation requirements without needing to be recompiled or modified as the scene changes. They can be adapted to different water wave simulation requirements simply by replacing the external perturbation input.

[0038] Step S312: Select the corresponding initial disturbance source data according to the type of the target water wave scenario. Among them, the earthquake tsunami scenario corresponds to the earthquake fault parameter data, the tidal bore scenario corresponds to the time-varying boundary data generated by the boundary condition function, and the landslide lake tsunami scenario corresponds to the landslide motion parameter data.

[0039] After the main simulation framework is solidified, the initial disturbance source data is selected according to the type of the target water wave scenario. This invention pre-defines the correspondence between three typical water wave scenarios and the initial disturbance source data. When an earthquake-tsunami scenario is designated as the simulation target, the initial disturbance source data is selected as earthquake fault parameter data. This data usually exists in a data file independent of the main program code, containing parameters describing the geometry and kinematic characteristics of the rupture surface, such as the epicenter location, fault strike, dip angle, slip angle, fault length, fault width, and slip amount. In practice, users only need to open the data file and modify the corresponding parameter fields line by line to complete the setting of the new earthquake scenario without having to enter the solver source code to make any changes. When a tidal bore scenario is designated, the initial disturbance source data is generated by the boundary condition function. This function is compiled into part of the executable program and provides time-varying water level or velocity values ​​to the grid nodes at the boundary of the computational domain at each time step to simulate the process of tidal waves entering the estuary or the propagation of tidal bores in the river channel. By changing the internal definition of the boundary condition function, it is possible to flexibly simulate boundary inputs driven by regular sinusoidal tidal waves, irregular composite tidal waves, or measured tidal process lines. When a landslide-tsunami scenario is specified, the initial disturbance source data corresponds to landslide motion parameter data. This data is jointly controlled by the landslide mass distribution function and landslide velocity function preset in the program, and is used to describe the mass transfer and kinetic energy input of the landslide body from initiation to sliding into the water body. The landslide effect is then transferred to the water body through the source term of the momentum equation. The correspondence between these three scenarios and data types is embedded in the loading process in the form of logical branches. The scenario identifier parameter controls the path selection, realizing branching operations under a unified entry point.

[0040] Step S313: Load the selected initial disturbance source data into the computation grid.

[0041] After completing scene determination and data type selection, the selected initial disturbance source data is loaded into the computational grid. For earthquake-tsunami scenarios, the loading process first reads the earthquake fault parameter data file, analyzes the rupture source location and geometric parameters of each fault marked in grid coordinates, then converts the fault slip amount into coseismic deformation distribution of the seabed based on the Okada fault model, and then converts the seabed deformation into the initial free water surface disturbance field according to the spatial correspondence between water surface and seabed deformation. Finally, the water surface disturbance value is assigned to each computational grid node near the rupture source, while the initial water level value at other locations in the computational domain is set to zero, and the initial velocity field value is set to zero across the entire domain, completing the mapping of the initial tsunami conditions onto the computational grid. For tidal bore scenarios, the loading process involves activating the boundary condition function call interface at the boundary of the computational domain. At the initial moment, a constant water level field and a zero velocity field can be preset throughout the entire computational domain. The boundary condition function gradually injects water movement into the computational domain over time, and the time-varying water level or velocity values ​​on the boundary grid are synchronized to the boundary nodes of the computational grid according to the time step. For the landslide-lake tsunami scenario, a kinematic landslide model is used for simulation. During the loading process, the landslide mass distribution function and velocity function are associated with the grid of the specified landslide area within the computational domain. Initially, the water is stationary. As time progresses, the landslide kinematic function outputs the sliding velocity and mass flux of the landslide at different spatial locations step by step, and these are loaded as the volumetric force source terms of the momentum equation into the discrete equations of the corresponding grid nodes. The loading operations for all three scenarios do not modify the main solver structure within the hybrid numerical discretization framework; scenario adaptation is only completed at the level of external data mapping and source term assignment.

[0042] Step S314: When the target water wave scene changes, replace the type of the initial disturbance source data and the corresponding source item processing operation, and return to execute step S312.

[0043] When a simulation task needs to switch from one water wave scenario to another, the scenario switching process is triggered by resetting the scenario identifier parameters. This process first clears the initial disturbance field and boundary condition settings left over from the previous scenario, including resetting the water surface disturbance field, releasing boundary condition function pointers, or clearing intermediate state variables of the landslide motion function. Then, based on the newly specified scenario identifier, it replaces the type of the initial disturbance source data and the corresponding source item processing operations: if the new scenario is an earthquake-tsunami, it re-specifies the earthquake fault parameter data file path and calls the fault-to-water surface disturbance mapping program; if the new scenario is a tsunami, it recompiles and links the corresponding boundary condition functions; if the new scenario is a landslide-lake tsunami, it re-specifies the landslide mass distribution function and motion velocity file. After completing the data type replacement, the program returns to the scenario determination and data selection steps, and re-loads the initial disturbance source data and establishes the initial conditions on the computational grid according to the new scenario configuration. Throughout the scene switching process, the main solver structure of the simulation framework remains unchanged. The mesh topology, discretization format, time progression logic, and wet / dry determination rules do not need to be recompiled or modified. Only the input level of the disturbance source is unloaded and remounted. This mechanism enables three long-wave disaster problems with completely different causes—tsunamis, tidal bores, and landslide tsunamis—to be switched and analyzed continuously on the same numerical solution platform as needed.

[0044] By constructing a simulation framework and solidifying its solver structure, a standardized correspondence between three typical water wave scenarios and initial disturbance source data is established. The data selection and loading process is driven by scenario identifier parameters. At the same time, the solver structure remains unchanged when switching scenarios, and only the disturbance source type and source term processing operation are replaced. This allows three long-wave problems with different causes—tsunamis, tidal bores, and landslide tsunamis—to be numerically simulated within a unified hybrid numerical discretization framework. This eliminates the code redundancy and maintenance costs associated with maintaining dedicated solvers for different scenarios. It also provides a technically consistent numerical experimental platform for the systematic comparison of water wave propagation and evolution processes under different disaster-causing mechanisms. The overall solution has good engineering versatility and scenario expansion capabilities.

[0045] Furthermore, the initial disturbance source data is seismic fault parameter data, which is stored in a data file format independent of the solver code.

[0046] When the initial disturbance source data is concretized into seismic fault parameter data, this data is stored and managed as a data file independent of the solver code. This means that all parameters describing the characteristics of earthquake rupture (including focal location, fault strike, dip angle, slip angle, fault length, fault width, and slip amount, etc.) are not hard-coded into the solver's main program source code, but are stored in a separate editable text file. When simulating different seismic events, users only need to open this data file with a text editor, modify the corresponding parameter values ​​line by line or adjust the wave source location field marked in grid coordinates, save, and then rerun the solver to complete the setting of the new scenario, without needing to make any changes or recompile at the source code level. This engineering implementation method of separating physical parameters from the computational kernel lowers the operating threshold for non-programming users, avoids program logic errors that may be introduced by directly modifying the code, and allows the same solver executable program to quickly adapt to tsunami scenarios under different focal mechanisms by simply replacing the data file, significantly improving the reusability of simulation tools in multi-scenario batch calculations and parameter sensitivity analysis.

[0047] Accordingly, step S313, loading the initial disturbance source data into the computational grid corresponding to the governing equations, includes: Step S3131: Read the data file, obtain the wave source position parameters specified in the data file in the form of grid coordinates, and map the wave source position to the computation grid.

[0048] During the loading of earthquake fault parameter data into the computational grid, the program first reads a data file stored independently of the solver code. This data file is in plain text format, storing parameters such as the epicenter location, fault strike, dip angle, slip angle, fault length, fault width, and slip amount line by line. The epicenter location is directly labeled in grid coordinates; that is, the user writes the row and column indices of the target rupture source in the current computational grid, rather than geographic latitude and longitude. This grid coordinate labeling of the epicenter location is designed to directly interface with the solver's internal spatial data structure, which uses grid cells as the basic index unit, avoiding repeated conversions between geographic and grid coordinates during the loading process. During reading, the program parses the file content line by line, storing each parameter field into the corresponding variable structure. Then, it extracts the grid row and column indices, locates the corresponding grid cell within the two-dimensional terrain array and the initial water surface disturbance array in memory, and delineates the spatial distribution area of ​​the fault rupture surface centered on this cell. In practical implementation, for example in a global tsunami simulation scenario, if the current computational grid size is 1440×720, and the grid coordinates that can be used to mark the location of the rupture source in the data file are the 720th grid in the meridional direction and the 360th grid in the zonal direction, after the program reads the data, it establishes the initial disturbance center at the grid location and maps the dislocation amounts on both sides of the rupture surface to the initial water level values ​​of the surrounding nodes according to the fault strike and dip parameters, thus completing the spatial mapping of the wave source location from the file data to the computational grid.

[0049] Step S3132: When the resolution of the computation grid changes, the geographical location of the wave source in the target area is kept unchanged by adjusting the grid coordinate parameters in the data file.

[0050] When the resolution of the computational grid is adjusted due to changes in simulation accuracy requirements, the wave source location parameters specified in grid coordinates in the data file need to be corrected accordingly to ensure that the actual geographical location of the wave source in the target area does not shift due to grid densification or sparseness. This resolution adjustment is quite common in practical applications: for example, in the initial calculation phase of a global tsunami simulation, a coarser resolution grid of 15 arcminutes might be used to quickly obtain the basic characteristics of wave propagation over a large area; however, when performing refined wavefield analysis on a specific nearshore area, the global grid resolution needs to be densified to 10 arcminutes or even higher. In this case, the rupture source under the same geographical coordinates will fall into different grid index numbers in different grids. If the grid coordinate values ​​labeled for the coarse grid in the original data file are continued to be used, the geographical location of the wave source will shift westward or southward in the new grid, causing a systematic deviation between the simulated initial tsunami disturbance field and the actual earthquake epicenter location. To solve this problem, this step allows users to complete the correction by editing the wave source grid coordinate field in a separate data file without modifying the solver source code. The specific steps are as follows: Recalculate the row and column numbers of the rupture source location in the new grid based on the ratio of the old to the new resolution, and update the corresponding field values ​​in the data file. For example, when the meridional grid resolution is increased from 15 arcminutes to 10 arcminutes (meaning the span of each grid cell is reduced to 2 / 3 of its original size, and the number of grid cells increases to 1.5 times), multiplying the meridional grid coordinate values ​​in the data file by a scaling factor of 1.5 will ensure that the wave source location falls on the same geographic coordinate point in the new grid as before. After adjusting the grid coordinate parameters, save the data file and rerun the solver to obtain simulation results that are geographically consistent with the resolution change.

[0051] By storing earthquake fault parameters as independent data files and directly labeling the wave source location in the files using grid coordinates, and by supporting the modification of grid coordinate values ​​in the files to keep the geographical location of the wave source unchanged when the grid resolution changes, the same solver executable program and the same parameter management system can stably adapt to different computational needs from coarse-resolution rapid assessment to high-resolution fine simulation. Both the engineering level of parameter management and the geographical consistency of simulation results are guaranteed.

[0052] Furthermore, the initial disturbance source data consists of earthquake fault parameter data, including: focal location, fault strike, dip angle, slip angle, fault length, fault width, and / or slip amount. The focal location is marked with grid coordinates or geographic coordinates indicating the spatial position of the rupture initiation point in the computational domain. The fault strike and dip angle together describe the attitude of the fault plane in three-dimensional space. The slip angle defines the slip direction of the hanging wall relative to the footwall. The fault length and width determine the spatial distribution of the rupture surface, and the slip amount quantifies the relative dislocation magnitude between the two sides of the fault. This set of parameters constitutes a complete parameter set describing the geometric and kinematic characteristics of the earthquake rupture. After being read by the program, the parameters are converted into coseismic deformation distributions of the seabed surface according to elastic dislocation theory, and then mapped to the initial water surface disturbance field of the governing equations, providing initial conditions consistent with the actual earthquake mechanism for the numerical simulation of the tsunami propagation process.

[0053] Specifically, in step S3, the conservation terms are discretized using the finite volume method, including: Step S321: Extract the physical quantities of the mesh cells on both sides of the interface and calculate the local gradient ratio of the physical quantities.

[0054] When discretizing conservation terms using the finite volume method, the physical quantity information of adjacent grid cells on both sides of each grid interface needs to be extracted first, and the local gradient ratio is calculated accordingly. This physical quantity is typically a water level or velocity component, depending on the type of conservation equation being solved. After extracting the physical quantities on both sides of the interface, the projection of their spatial gradient onto the interface normal direction is calculated, and the gradient on one side of the interface is compared with the gradient on the other side to obtain a dimensionless local gradient ratio parameter. This parameter essentially reflects the intensity of flow field changes near the interface: when the local gradient ratio is close to 1, it indicates that the physical quantities on both sides of the interface change gradually, and the flow field is in a smooth region; when the local gradient ratio deviates significantly from 1, it indicates strong gradient changes or even discontinuous structures near the interface, such as the abrupt change in water level within several grid intervals when a tidal wave front advances in a river channel. By calculating this local gradient ratio on each grid interface, a quantitative local smoothness criterion can be provided for the subsequent adaptive adjustment of the limiter function, enabling the numerical scheme to adopt different interpolation strategies in smooth and discontinuous regions.

[0055] Step S322: Input the local gradient ratio into the flux limiter function, and the flux limiter function determines the interface interpolation weights used to calculate the flux.

[0056] The calculated local gradient ratio is passed as an input parameter to the flux limiter function, which dynamically determines the interface interpolation weights used to calculate the interface flux based on the magnitude of the local gradient ratio. The flux limiter function is a predefined nonlinear mapping relationship; its input is the local gradient ratio, and its output is the interpolation weight coefficients. When the local gradient ratio indicates that the interface is in a smooth flow field region, the interpolation weights output by the limiter function tend to use a higher-order interpolation scheme, such as assigning weights with larger central differences, to take advantage of the high accuracy of higher-order schemes in smooth regions. When the local gradient ratio indicates the presence of strong gradients or discontinuous structures near the interface, the interpolation weights output by the limiter function automatically tilt towards a lower-order upwind scheme, sacrificing some accuracy for numerical stability. In practical implementation, the flux limiter function can support various preset mathematical forms, and users can select the appropriate limiter type before calculation based on the gradient characteristics of the simulated scenario. This mechanism, which adaptively adjusts the interpolation weights based on the actual local flow field conditions, eliminates the need for manual judgment of smoothness at each interface. Instead, the limiter function performs the calculation uniformly, enabling the entire flux reconstruction process to automatically adapt to the discrete needs of different regions in a complex flow field within a unified mathematical framework.

[0057] Step S323: Based on the interface interpolation weight, the physical quantities on both sides of the interface are reconstructed by weighting to obtain the interface reconstruction value used to calculate the flux.

[0058] Based on the interface interpolation weights determined in the previous step, the physical quantities of the grid cells on both sides of the interface are reconstructed using weighted averages to obtain the interface reconstruction values ​​used in the flux calculation formula. The weighted reconstruction operation uses the physical quantities at the center of the grid on both sides of the interface as a basis, and linearly combines them according to the interpolation weights. This extrapolates or interpolates the physical quantities originally defined at the grid center point to the center point of the grid interface, forming a one-sided reconstruction value or a two-sided average value on the interface. In smooth flow fields, because the interpolation weights are close to the central difference mode, the reconstruction values ​​can retain the higher-order spatial variation information of the physical quantities near the interface, making the flux calculation more accurately reflect the actual physical transport. In strong gradient regions, the interpolation weights are biased towards the windward direction, and the reconstruction values ​​are closer to the physical quantities at the center of the upwind grid, thereby suppressing non-physical numerical overshoot caused by sharp jumps between the two sides of the interface. For example, at the grid interface through which a tidal wavefront passes, the water level on either side of the interface may differ by several meters. If central difference reconstruction is used indiscriminately, severe numerical oscillations will occur. However, the interpolation weights adjusted by the limiter function will cause the reconstructed values ​​to preferentially accept the water level information upstream of the wavefront, ensuring that flux calculations remain bounded even under extreme hydraulic gradient conditions. The interface physical quantities after weighted reconstruction are the interface reconstruction values, which serve as the direct input for subsequent flux calculations.

[0059] Step S324: Substitute the interface reconstruction value into the flux calculation formula to calculate the flux through each grid interface, and use the flux to update the conserved quantity.

[0060] Substituting the interface reconstruction values ​​obtained through weighted reconstruction using limiters into the corresponding flux calculation formulas, the numerical fluxes through each grid interface are calculated. The choice of flux calculation formula depends on whether the current solution is a continuity equation or a momentum equation: for the continuity equation, the flux is usually the product of the interface reconstruction velocity and the interface area; for the momentum equation, the flux includes contributions from convective transport and pressure terms. The interface flux obtained by substituting the interface reconstruction values ​​into the formulas represents the mass or momentum transport through that interface per unit time. After completing the flux calculations for all grid interfaces, a conservation update is performed on each grid cell, i.e., the interface fluxes flowing into and out of the cell are algebraically summed to obtain the net flux increment of the cell in the current time step, and then the conservation variables of the cell, such as water level or momentum components, are updated. Through this conservation update method based on interface reconstruction fluxes, the governing equations strictly satisfy the physical conservation laws at the discrete level. The changes in physical quantities within each grid cell are entirely determined by the net flux input through its boundaries, ensuring the numerical conservation of total mass and total momentum across the entire domain.

[0061] By introducing a flux limiter function based on the local gradient ratio into the discretization of the conservation terms, the physical quantities at the interface are adaptively weighted and reconstructed. The reconstructed values ​​are then substituted into the flux formula to update the conservation quantities. This allows the numerical scheme to automatically approach high-order accuracy in smooth flow fields to accurately characterize the details of water wave propagation, while in strong gradient and discontinuous regions it automatically reduces to a low-order upwind mode to effectively suppress non-physical numerical oscillations. Thus, within a unified finite volume method framework, both high accuracy and high stability of flux calculation are achieved.

[0062] Furthermore, the flux limiter function can be upwind, QUICK, or Superbee format. The flux limiter function uses higher-order interpolation weights in low-gradient regions and lower-order interpolation weights in high-gradient regions, based on the magnitude of the local gradient ratio.

[0063] In the specific implementation of the flux limiter function, its selectable format is further limited to one of the upwind, QUICK, or Superbee formats. The upwind format only uses the physical quantities of the upwind grid when calculating the interface interpolation weights, possessing first-order accuracy and natural monotonicity; the QUICK format introduces second-order curvature corrections of the physical quantities on both sides of the interface along the flow direction on the upwind basis, achieving third-order accuracy in smooth regions; the Superbee format is a high-resolution limiter, striking a better balance between intermittent sharpness capture and suppression of numerical diffusion. Regardless of the format chosen, the core adjustment mechanism is consistent: the flux limiter function uses the local gradient ratio as a criterion, automatically increasing the weights from the central difference or higher-order template in the interface interpolation in low-gradient regions with gentle flow field changes, in order to approximate higher-order accuracy in these regions; in high-gradient regions caused by wavefronts or hydraulic abrupt changes, it automatically reduces the weights of higher-order contributions and shifts to lower-order interpolation weights that are biased towards the upwind side, thereby effectively suppressing numerical overshoot and spurious oscillations at these locations.

[0064] Specifically, in step S3, the dispersion term is discretized using the finite difference method, including: Step S331: Split the spatiotemporal mixed partial derivatives contained in the dispersion term, transform the mixed term of the second-order time derivative and spatial derivative into an algebraic combination of the first-order spatial derivative and the second-order spatial derivative, and use a local semi-implicit method to advance the time.

[0065] When discretizing dispersion terms using the finite difference method, the spatiotemporal mixed partial derivatives contained in the dispersion terms are first split into terms. Mathematically, dispersion terms in Boussinesq-type equations are usually represented as mixed higher-order partial derivative terms containing the product or coupling of the second-order time derivative and the spatial derivative, such as the combination of the second-order time derivative of water level with respect to space. These mixed derivatives cannot be directly discretized using the standard difference template. The core idea of ​​the spatiotemporal hybrid partial derivative term decomposition is to extract the contribution of the central node related to the momentum variable to be solved in the current grid from the dispersion term and add it to the left side of the update equation to form a local implicit stabilization coefficient. The dispersion time increment of adjacent nodes is explicitly approximated by variables of known time layer or lag time layer. In this way, the dispersion correction problem that originally required solving a global sparse linear equation system is transformed into a local algebraic update process per grid point. At the same time, by combining the finite volume method to discretize the conservation flux term, the finite difference method to discretize the higher-order dispersion term, and dynamic CFL time step control, the computational complexity and memory overhead under the condition of large-scale spherical grid are significantly reduced while ensuring the computational stability of flood inundation and tsunami propagation.

[0066] Step S332: On the computational grid, differentially discretize the higher-order spatial derivatives in the dispersion term and construct a local semi-implicit dispersion correction expression.

[0067] After transforming the dispersion term into a locally semi-implicit combination, the spatial first and second derivatives are numerically approximated on the constructed computational grid using the differences in physical quantities at adjacent nodes. The spatial first derivative employs a central difference scheme, approximating the partial derivative value at that point by the ratio of the difference in physical quantities at the center points of adjacent grids on both sides of the interface to the grid spacing. The spatial second derivative is approximated by further differencing the first derivative, or by directly approximating it with a weighted combination of physical quantities from three adjacent nodes, ensuring that the truncation error between the discrete result and the continuous partial derivative is controlled to the order of the square of the grid spacing. Through this difference approximation step, the dispersion term, originally represented by continuous partial derivatives, is transformed into a computable expression composed of the values ​​of water level or flow velocity components at adjacent grid nodes through addition, subtraction, multiplication, and division. Since all physical quantities in this expression come from the water level field and flow field known at the current time step, its calculation result can be directly used as an explicit contribution term in the numerical solution, without needing to incorporate the dispersion term into the coefficient matrix for simultaneous solution as in traditional methods. In the staggered grid configuration, the water level is defined at the center of the grid cell, and the velocity component is defined at the center of the grid boundary surface. The above differential operation can make full use of the natural alignment characteristics of each physical quantity at each custom position, avoiding additional interpolation propagation errors.

[0068] Step S333: Within the current time step, substitute the computable expression as a known source term into the governing equation to calculate the value of the dispersion term.

[0069] Within the current time step, the computable expression of the dispersion term obtained through differential approximation is directly substituted into the governing equations as a known source term for solution. Since the dispersion term, after decomposition and differential approximation, is fully expressed as the algebraic result of the physical quantities at the current moment, it no longer implicitly depends on unknown quantities at the next moment during time progression. Therefore, its value can be regarded as an external forcing term in the momentum equation or continuity equation, parallel to the pressure gradient term, convection term, etc., calculated once in the same explicit time progression step and superimposed on the intermediate state variables. In the program implementation, the calculation of the dispersion term is encapsulated in an independent sub-function, which is called in each time loop. The input is the current water level field and velocity field array, and the output is an array of discrete values ​​of the dispersion term at each grid node. This array is then directly accumulated to the total momentum change in the momentum update step. The entire solution process does not construct a system of matrix equations, does not execute iterative loops, and does not allocate additional storage arrays for the dispersion term beyond the lifetime of the current time step. This approach allows the solution of the most complex dispersion terms in high-order Boussinesq equations to be completed within a unified explicit framework, maintaining the simplicity and scalability of the overall program structure.

[0070] By splitting the dispersion term into terms using spatiotemporal mixed partial derivatives and employing local semi-implicit approximation, and then approximating the spatial derivatives of each order on the computational grid using the difference values ​​of adjacent nodes to form a computable expression, this expression is finally substituted into the governing equation as a known source term to directly calculate the value of the dispersion term within the current time step. This allows the calculation of higher-order dispersion contributions, which originally required implicit iteration or the introduction of auxiliary variables, to be completed within a unified explicit time-progression framework. This achieves an efficient numerical description of the dispersion effect of non-hydrostatic pressure without increasing program complexity or storage overhead.

[0071] Furthermore, the adjustable dispersion parameter participates in the construction of the computable expression in step S333 as a scaling factor. When the computable expression is substituted into the governing equation as a known source term in step S333, the adjustable dispersion parameter directly participates in the construction of this computable expression as a scaling factor. Specifically, in the algebraic combination expression obtained by term splitting and difference approximation of the dispersion term, the adjustable dispersion parameter is embedded as a dimensionless multiplier factor in the coefficients of each component of the dispersion term, uniformly scaling the numerical values ​​of the spatial derivative combination components of the dispersion term through multiplication. This parameter is not a physical constant, but an adjustable value pre-set by the user in the parameter configuration file according to the dispersion characteristics requirements of the simulation scenario before calculation. When the adjustable dispersion parameter takes different values, the values ​​of the dispersion term calculated under the same discretization scheme and the same set of flow field variables will be proportionally amplified or reduced, thereby directly adjusting the contribution ratio of the dispersion effect in the momentum equation. This method of embedding in the form of scaling factors is the simplest to implement. It does not require modifying the spatial derivative combination structure of the dispersion term, nor does it require introducing conditional branch statements in the program code. Continuous control of the dispersion intensity can be achieved simply by changing an external parameter.

[0072] When the adjustable dispersion parameter takes different values, the contribution intensity of the dispersion term in the governing equation is adjusted accordingly, allowing the governing equation to switch between the Peregrine dispersion form and the improved dispersion relation model. By assigning different specific values ​​to the adjustable dispersion parameter, the same governing equation framework can be smoothly switched between the classical Peregrine dispersion form and the improved dispersion relation model. When the adjustable dispersion parameter is zero, the scaling factor of the dispersion term is zero, the dispersion contribution in the equation disappears completely, and the governing equation degenerates into the basic Boussinesq-type equation dispersion form given in the classical derivation of Peregrine. When the adjustable dispersion parameter is a specific non-zero value, such as a rational number like 1 / 15, the dispersion term is activated proportionally and adjusted to a specific intensity. At this time, the linear dispersion relation curve of the equation will shift from the Peregrine form to a direction that more closely approximates the dispersion relation of the complete Airy wave theory, corresponding to the improved dispersion model optimized by the dispersion relation. This switching is achieved entirely by modifying a single value in the parameter configuration file. The dispersion term calculation sub-function in the solver's main code does not need to be recompiled or have its algorithm branches replaced as the parameters change. After changing the adjustable dispersion parameter from one value to another in the parameter file and saving it, the user can obtain simulation results under the new dispersion relationship by rerunning the executable program. This provides operational convenience for flexibly selecting the most suitable dispersion approximation accuracy under different water depth ranges and wavelength scales.

[0073] Specifically, step S4 involves truncating the dispersion term in the region involving dry mesh elements, including: Step S41: Before performing differential solution of dispersion term on the current wet mesh cell, determine all adjacent mesh cells covered by the dispersion term differential template.

[0074] In the computation of a hybrid numerical discretization framework, before the program performs a difference calculation of the dispersion term on a cell marked as a wet grid, it needs to first determine the range of all adjacent grid cells covered by the dispersion term difference template at the current grid location. After spatial discretization, the dispersion term is typically calculated using a multi-point centered difference template. Taking the difference of the spatial second derivative as an example, it requires physical quantity values ​​spanning at least three grid nodes in each coordinate direction to form a complete difference expression. If the dispersion term includes a spatial hybrid derivative, the difference template will also cover adjacent cells in the diagonal direction. Therefore, the grid cells relied upon for a dispersion term solution operation are not limited to the current computation location itself, but rather radiate outward from it to a template region of finite size. In actual execution, the program, based on the pre-set difference template radius, extracts all cell numbers from the grid array within a range of 1 to 2 grid indices extending outward in each coordinate direction from the current wet grid cell, forming a list of adjacent cells. For example, in a two-dimensional staggered grid, a standard second-order central difference template extends by one cell in both the longitudinal and latitudinal directions. The covered adjacent cells include the current cell itself and its four immediate neighbors to the east, south, west, and north. Determining this coverage area is a prerequisite for subsequent dry / wet state checks and truncation determination, ensuring that no adjacent cell that might cause numerical anomalies due to being located at the dry / wet boundary is overlooked.

[0075] Step S42: Check the dry / wet status indicators of all adjacent grid cells one by one.

[0076] After determining all adjacent grid cells covered by the differential template, the wet / dry status flags assigned to these cells in the previous time step or the current substep are checked one by one. The wet / dry status flag is an integer marker assigned after each time step water level update by comparing the current water depth of each grid cell with a preset water depth threshold; for example, 0 represents a wet cell and 1 represents a dry cell. Checking one by one involves reading the wet / dry flag of each cell sequentially along the adjacent cell list and determining whether it equals the flag value of a dry cell. In implementation, this detection process is compiled into a fast pre-loop at the entry point of the dispersion term solver function. As long as the number of adjacent cells is determined by the differential template radius, the number of loops is a fixed and finite constant, and no additional computational burden is introduced due to an increase in the total grid size. Taking the tidal wavefront climbing process as an example, when the wavefront reaches the gentle slope shoreline, the shore-side cells adjacent to the wet grid cells at the forefront of the wavefront are marked as dry grids because the water depth is zero. When the detection program traverses to these shore-side adjacent cells, it will identify the existence of dry grids, thereby triggering subsequent truncation logic. This unit-by-unit traversal detection mechanism ensures that no matter how complex the geometry of the wet-dry boundary line is in the mesh, the high-frequency discrete term solution region can be effectively identified at the entrance.

[0077] Step S43: If any dry grid cell exists among all adjacent grid cells, set the calculated value of the dispersion term at the current wet grid cell to zero and skip the subsequent solution steps for the dispersion term.

[0078] If any dry grid cell is found among all adjacent grid cells during the process of sequential detection, the dispersion term calculation at the current wet grid cell is immediately truncated. Specifically, the truncation involves forcibly setting the dispersion term calculation value of that wet grid cell to zero at the current time step and skipping subsequent solution steps for the dispersion term at that cell, including skipping the differential summation of the spatial derivative and the subsequent process of accumulating the result to the source term of the momentum equation. This processing is based on the following considerations: the dispersion term physically describes the nonlocal contribution of the non-hydrostatic pressure effect of water waves to momentum. When there are adjacent dry grid cells in the differential template, these cells lack a physical definition of water, and the corresponding physical quantities are zero or undefined. Forcibly including them in the differential summation would mathematically introduce a false one-sided derivative jump, equivalent to artificially injecting non-physical momentum perturbations at the water-land transient interface. Setting the dispersion term to zero is a conservative but numerically safe approach. Physically, it means ignoring the dispersion effect of fluctuations at the wet-dry interface and retaining only the convection and pressure contributions of the shallow water equations. This aligns with the physical tendency for the dispersion effect in shallow water to naturally weaken due to the extremely shallow water depth. After truncation, the program will mark the skipped state and move on to the next wet grid cell, without returning to perform dispersion calculations on the skipped cells.

[0079] Step S44: If all adjacent mesh elements are wet mesh elements, then calculate according to the differential combination expression of the dispersion term.

[0080] If, after individual checks, all adjacent grid cells are confirmed to be wet grid cells (meaning the water depth of each grid cell covered by the difference template is greater than a preset threshold and the wet / dry status identifiers are all wet cell markers), then the current wet grid cell is located in the inner water area, and there are no boundary cells with missing data around it. At this point, the program follows the normal calculation path, extracting the required physical quantities from each adjacent cell according to the difference combination expression determined after term splitting of the dispersion term. It then calculates the difference approximations of the first and second spatial derivatives, substitutes them into the computable expression to obtain the discrete result of the dispersion term in that cell, and normally superimposes this result into the dispersion term contribution array of the momentum equation. In the inner water area, the water depth of all adjacent cells has physical meaning, the difference operation is isotropic, and the numerical approximation accuracy of the dispersion term is fully preserved. This normal calculation path and the aforementioned truncated path together constitute a two-branch logical structure within the dispersion term solution subfunction. The wet / dry detection results of adjacent cells serve as the branch judgment condition, automatically selecting the execution path for each wet grid cell at each time step, without user intervention or manual division of the computational domain for different terrain regions.

[0081] By determining the range of all adjacent cells covered by the differential template before solving the dispersion term and checking the dry / wet status of each cell one by one, the dispersion term is forcibly set to zero and the solution is skipped when any adjacent cell is a dry grid. The normal difference expression is only used when all adjacent cells are wet grids. This ensures that non-physical numerical disturbances that may be caused by missing data in the dry / wet boundary region are truncated and isolated at the entrance. The numerical stability of the dispersion term solution is systematically guaranteed in the entire computational domain. The reliability and physical rationality of the simulation results under complex shoreline dynamic boundary conditions are significantly improved.

[0082] like Figures 2-9 As shown, Figures 2-9 This simulation encompasses the entire process from calm global ocean surfaces to the generation and global propagation of tsunamis triggered by earthquakes. It demonstrates the generation and propagation evolution of tsunamis in global ocean areas, intuitively reflecting the spatial response relationships of wave signals at different times between land boundaries, island arcs, ocean basins, and nearshore regions. It shows that the model can uniformly visualize the calculation results of global-scale topography and tsunami wave fields, thus providing intuitive evidence for determining the impact range of transoceanic tsunamis, identifying key coastal risks, analyzing arrival time, disaster warning, and emergency response.

[0083] like Figures 10-13 As shown, Figures 10-13 It shows the initial state of the Qiantang River, the generation and propagation of the tidal bore, and intuitively demonstrates the spatial distribution of water surface elevation at different times, the range of tidal wave propagation, the changes in water level near the shore and in the river channel, and the relationship between the humid and dry areas. It can intuitively illustrate that the model can simulate the propagation, deformation, slope climbing and inundation capabilities of tidal waves in complex estuary topography.

[0084] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor and a memory connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, cause the at least one processor to perform the aforementioned numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations.

[0085] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations.

[0086] The embodiments of this invention aim to protect a numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations, which has the following effects: 1. By splitting the spatiotemporal mixed partial derivatives in the dispersion term and transforming them into an algebraic combination of the first and second spatial derivatives, this invention constructs a locally semi-implicit time-propelling architecture. This allows higher-order dispersion terms to be directly substituted into the governing equations as known source terms. At the same time, adjustable dispersion parameters are introduced as scaling factors to participate in the calculation of the dispersion term. Thus, under a unified explicit time-propelling framework, the approximation accuracy of the dispersion relationship can be dynamically adjusted without iterative solutions. This achieves an accurate description of the propagation characteristics of different wavelength components while significantly reducing computational complexity. 2. By introducing a dynamic identification mechanism for wet and dry grids based on a preset water depth threshold into the hybrid numerical discretization framework, and by checking the state of adjacent grid cells covered by the differential template one by one before calculating the dispersion term, the solution of the dispersion term is truncated once a dry grid cell is found. This effectively curbs the non-physical numerical oscillations caused by data loss in the wet and dry interface area, and ensures the stability and reliability of the simulation results of the flow field and water level field in the whole computational domain during the advance and fall of the water body. 3. By constructing a simulation framework composed of governing equations, a hybrid numerical discretization framework, and dynamic processing of wet and dry grids, and storing earthquake fault parameters in independent data files or defining initial disturbance sources for different scenarios using boundary condition functions and landslide motion parameters, the same solver structure can maintain the geographical location of the wave source at different grid resolutions and flexibly switch between different water wave scenarios. Combined with the flux limiter reconstruction of interface flux in the finite volume method discretization, numerical overshoot and dissipation in high gradient regions are further suppressed. Overall, high-precision and robust engineering simulation of various complex hydrodynamic processes such as tsunamis, tidal bores, and landslide tsunamis is achieved.

[0087] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0088] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0089] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0090] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0091] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations, characterized in that, Includes the following steps: Step S1: Obtain digital elevation data and initial disturbance source data for the target area; Step S2: Construct a control equation containing an adjustable dispersion parameter. The control equation is a high-order Boussinesq type equation that introduces a dispersion term. The adjustable dispersion parameter is used to adjust the approximation accuracy of the control equation for the difference in propagation speed of different wavelength components. Step S3: Establish a hybrid numerical discretization framework, load the digital elevation data and the initial disturbance source data into the computational grid corresponding to the control equation, and discretize the conservation terms in the control equation using the finite volume method and the dispersion terms in the control equation using the finite difference method under a unified explicit time-progression system. Step S4: During the calculation process of the hybrid numerical discretization framework, a dry-wet grid dynamic processing is performed to obtain the numerical simulation results of the flow field and water level field. The dry-wet grid dynamic processing includes: determining the dry-wet state of the grid cell based on a preset water depth threshold, and truncating the solution of the dispersion term in the region involving the dry grid cell.

2. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 1, characterized in that, Step S3 involves discretizing the dispersion term using the finite difference method, including: Step S331: Split the spatiotemporal mixed partial derivatives contained in the dispersion term, and transform the mixed term of the second-order time derivative and spatial derivative into an algebraic combination of the first-order spatial derivative and the second-order spatial derivative. Step S332: On the computational grid, differentially discretize the higher-order spatial derivatives in the dispersion term and construct a local semi-implicit dispersion correction expression; Step S333: Within the current time step, substitute the computable expression as a known source term into the governing equation to calculate the value of the dispersion term.

3. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 2, characterized in that, The adjustable dispersion parameter participates in the construction of the computable expression in step S333 in the form of a scaling factor; When the adjustable dispersion parameter takes different values, the contribution strength of the dispersion term in the control equation is adjusted accordingly, so that the control equation switches between the Peregrine dispersion form and the improved dispersion relationship model.

4. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 1, characterized in that, Step S3 involves loading the initial disturbance source data into the computational grid corresponding to the governing equations, including: Step S311: The simulation framework is composed of the governing equations, the hybrid numerical discretization framework, and the dry and wet mesh dynamic processing. The solver structure of the simulation framework remains unchanged in the following steps. Step S312: Select the corresponding initial disturbance source data according to the type of the target water wave scene; wherein, the earthquake tsunami scene corresponds to the earthquake fault parameter data, the tidal bore scene corresponds to the time-varying boundary data generated by the boundary condition function, and the landslide lake tsunami scene corresponds to the landslide motion parameter data. Step S313: Load the selected initial disturbance source data into the computational grid; Step S314: When the target water wave scene is switched, replace the type of the initial disturbance source data and the corresponding source item processing operation, and return to execute step S312.

5. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 4, characterized in that, The initial disturbance source data is seismic fault parameter data, which is stored in a data file format independent of the solver code; Step S313, loading the initial disturbance source data into the computational grid corresponding to the control equation, includes: Step S3131: Read the data file, obtain the wave source position parameters specified in the data file in the form of grid coordinates, and map the wave source position to the computation grid; Step S3132: When the resolution of the computational grid changes, the geographical location of the wave source in the target area is kept unchanged by adjusting the grid coordinate parameters in the data file.

6. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 5, characterized in that, The initial disturbance source data is earthquake fault parameter data, which includes: source location, fault strike, dip angle, slip angle, fault length, fault width, and / or slip amount.

7. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 1, characterized in that, Step S3 involves discretizing the conserved terms using the finite volume method, including: Step S321: Extract the physical quantities of the mesh cells on both sides of the interface, and calculate the local gradient ratio of the physical quantities; Step S322: Input the local gradient ratio into the flux limiter function, and the flux limiter function determines the interface interpolation weights used to calculate the flux; Step S323: Based on the interface interpolation weight, the physical quantities on both sides of the interface are reconstructed by weighting to obtain the interface reconstruction value used to calculate the flux; Step S324: Substitute the interface reconstruction value into the flux calculation formula to calculate the flux through each grid interface, and use the flux to update the conserved quantity.

8. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to claim 7, characterized in that, The flux limiter function is in upwind, QUICK, or Superbee format. The flux limiter function uses higher-order interpolation weights in low-gradient regions and lower-order interpolation weights in high-gradient regions, based on the magnitude of the local gradient ratio.

9. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to any one of claims 1-8, characterized in that, Step S1 involves obtaining the digital elevation data, including: Step S111: Obtain digital elevation data of the study area in its original format; Step S112: Convert the original format digital elevation data into a three-column data file of a unified format, wherein the three-column data file sequentially includes a longitude data column, a latitude data column, and a topographic elevation data column; Step S113: Using the three-column data file as the input format for the digital elevation data, load it into the computational grid.

10. The numerical simulation method for tsunamis and lake tsunamis based on Boussinesq-type equations according to any one of claims 1-8, characterized in that, Step S4, which involves truncating the solution of the dispersion term in the region involving dry mesh cells, includes: Step S41: Before performing the differential solution of the dispersion term on the current wet mesh cell, determine all adjacent mesh cells covered by the dispersion term differential template; Step S42: Detect the dry / wet status indicators of all adjacent grid cells one by one; Step S43: If any dry grid cell exists among all adjacent grid cells, set the calculated value of the dispersion term at the current wet grid cell to zero and skip the subsequent solution steps for the dispersion term. Step S44: If all adjacent mesh cells are wet mesh cells, then calculate according to the differential combination expression of the dispersion term.

Images