A method for constructing a three-dimensional spherical dynamic model of suboceanic plateau subduction

Abstract

The application discloses a kind of submarine plateau subduction three-dimensional spherical dynamics model construction methods, steps are as follows: solving mass, momentum, energy and composition conservation equation under the influence of neglecting compressible fluid dynamics;Apply mantle convection activity;Apply viscosity jump;Shape independent plate tectonics behavior;Embedding with multi-parameter attribute submarine plateau;After model test, finally complete the construction of three-dimensional spherical earth dynamics model.The three-dimensional numerical model constructed adopts three-dimensional spherical tectonic method with independent plate tectonics behavior, by introducing different sizes and different rheological state parameters, fine description is carried out to submarine plateau.The model can not only accurately reproduce the bending and fracture of plate and other geodynamic phenomena in the process of submarine plateau subduction, but also can reveal the difference of local trench tectonic response of different size and rheological property submarine plateau, and realize the cross-scale simulation from local structure to global structure influence.

Description

Technical Field

[0001] This invention relates to the technical fields of geophysics, geodynamics, marine geology, numerical simulation and computational geophysics, structural geology, and marine science. Background Technology

[0002] Given the difficulty and high cost of directly obtaining information about the deep Earth, numerical simulation has long been the mainstream international method for analyzing the dynamic mechanisms of submarine plateau subduction. Such simulations not only reproduce the dynamic evolution and deformation characteristics revealed by observational and survey data, but also overcome the limitations of geological time and space scales, enabling comprehensive discussions of various Earth science issues and providing quantitative support for observational results. Therefore, they have become an indispensable tool in marine geoscience research. Based on spatial dimensions, numerical simulations are generally divided into two categories: two-dimensional simulations and three-dimensional simulations.

[0003] In two-dimensional models, numerous studies have focused on the subduction process of oceanic plates carrying submarine plateaus. Yan et al. (2021) found through two-dimensional thermodynamic simulations that when a submarine plateau has a lithospheric basement of a certain thickness and the continental plate is undergoing thrusting, the insufficient buoyancy of the plateau can hinder subduction and ultimately trigger a leap. Dai et al. (2020)'s two-dimensional thermodynamic model also showed that large submarine plateaus on the Paleo-Pacific Plate could induce plate subduction, evolving into steep subduction from the Mesozoic to the Cenozoic. However, because such simulations assume no change in the direction of the trench, they highly simplify the geological processes, limiting their applicability. In contrast, three-dimensional numerical simulations, with their higher complexity and precision, can produce more realistic results closer to nature. For example, Liu et al. (2010), combining three-dimensional simulations with paleoplate reconstruction, pointed out that the timing of the subduction and removal of the Shatsky Submarine Plateau coincided with the Larami Orogeny, indicating that plateau subduction may trigger regional surface rebound and uplift. Betts et al. (2015) introduced the mantle plume effect into a three-dimensional experiment to further reveal the key influence of the synergistic effect of submarine plateaus and mantle plumes on the geometric evolution of trenches and subduction slabs.

[0004] However, current two-dimensional and three-dimensional studies on the subduction effect of submarine plateaus largely focus on ocean-continent convergence zones, with insufficient attention paid to processes within ocean-ocean convergence systems. Furthermore, existing three-dimensional models often neglect the influence of overlying plates and struggle to simultaneously integrate key parameters such as plate strength, mantle flow and thermal structure, and multi-rheological submarine plateaus. In particular, the relationship between the subduction of multi-sized, multi-rheological submarine plateaus and subduction zones has not been systematically explored. Therefore, current three-dimensional simulations still cannot fully represent the overall impact of submarine plateau subduction on the tectonic and dynamic behavior of subduction zones. Summary of the Invention

[0005] This invention discloses a method for constructing a three-dimensional spherical dynamic model of submarine plateau subduction, aiming to solve the problem that existing three-dimensional model studies neglect ocean-ocean plate convergence and overlying plate interactions. By using viscosity jump and pseudoplastic rheological approximations, a self-consistent plate-like tectonic behavior is generated, and a submarine plateau with multiple sizes and rheological parameters is introduced to realize the construction of a three-dimensional spherical dynamic model of submarine plateau subduction.

[0006] This invention is achieved through the following technical solution:

[0007] A method for constructing a three-dimensional spherical dynamic model of a subducting submarine plateau includes the following steps:

[0008] 1) Solve the conservation equations for mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids;

[0009] First, the high-viscosity flow equations are solved in three-dimensional space. Specifically, the governing equations are spatially discretized on an interlaced grid, ensuring that the different unknowns to be solved are not located in the same place. Pressure and temperature are located at the center of each element volume, while velocity components are located at the center of the corresponding element surface, avoiding artificial oscillations in the pressure field. The solver can handle relatively strong viscosity changes over small distances. To make the obtained model more realistic, viscosity can be an arbitrary function of temperature, depth, composition, and strain rate.

[0010] 2) Apply mantle convection activity

[0011] In the numerical model, the Rayleigh number is introduced as a measure of mantle convection activity;

[0012] 3) Apply viscosity jump

[0013] Viscosity can vary with temperature, depth, stress or strain rate, composition, phase transition, and melt fraction, increasing the viscosity at a depth of 660 km by 30 times;

[0014] 4) Shaping the behavior of plate-like structures

[0015] The model employs a pseudoplastic rheological approximation, enabling it to autonomously generate plate-like tectonic behavior, obtain Earth's surface velocity and tectonic activity, and form plate boundaries within and around rigid plates.

[0016] 5) Embedding submarine plateaus with multi-parameter attributes

[0017] Based on this, a global three-dimensional spherical mantle convection model was constructed. On this model, an ocean-to-ocean subduction zone was selected as the test area. A submarine plateau with multiple parameters was set up on the seaward side of the subducting ocean plate and close to the trench. This submarine plateau is characterized by multiple parameters such as geometric dimensions (length, width, and thickness) and rheological properties (buoyancy, viscosity, and yield stress).

[0018] 6) Test the model to verify the accuracy of the model solver in terms of temperature field, velocity field, buoyancy drive, etc., and finally construct a three-dimensional spherical seabed plateau subduction dynamic model.

[0019] Step 1) describes solving the conservation equations for mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids:

[0020] The variable viscosity is emphasized on the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude latitude and longitude grids. The mass conservation equation, momentum conservation equation, energy conservation equation and composition conservation equation in classical fluid dynamics are solved under the condition of considering the Businesk approximation, that is, neglecting the dynamic effects of compressible fluids.

[0021] Momentum conservation equation:

[0022]

[0023] mass conservation equation:

[0024]

[0025] Energy conservation equation:

[0026]

[0027] Composition conservation equation:

[0028]

[0029] in, Let Y represent the deviatoric stress tensor, Y represent pressure, R represent the Rayleigh number (a control parameter), b represent the radius, m represent density, F represent composition, and W represent total temperature. s B represents speed. r Ha represents specific heat capacity. b Surface dissipation number, r p The coefficient of thermal expansion, s b Represents radial velocity, r c N represents the thermal conductivity coefficient, and N represents the internal heating rate. Denotes the strain rate tensor, Δm the This represents the fractional density as a function of temperature.

[0030] Step 2) describes the application of mantle convection activity as follows:

[0031] The Rayleigh number is used to characterize the convective activity of the mantle in a model, and is defined as follows:

[0032]

[0033] Where G represents gravitational acceleration, m0 represents the reference density, and r p0 The coefficient of thermal expansion of the surface is represented by ΔW, where ΔW represents the temperature difference of the mantle region of thickness M, M represents the thickness of the mantle region, and r is the temperature coefficient of thermal expansion of the surface. k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

[0034] Step 3) describes the application of viscosity jump as follows:

[0035] 3.1) Baseline viscosity: Viscosity varies radially and laterally, increases exponentially with depth, and is also temperature-dependent, defined as follows:

[0036]

[0037] Where d represents depth, Wj represents absolute temperature, and A e The activation energy is represented by Q, and the gas constant is represented by Q.

[0038] 3.2) Viscosity jump: At the 660km phase boundary, the baseline viscosity is multiplied by the jump factor f. jump (e.g., 30×) to reflect the viscosity discontinuity caused by rock phase transformation. The viscosity with jump correction is then obtained:

[0039] Nd jump (d,Wj)=Nd(d,Wj)f jump (7).

[0040] The shaping behavior of the plate described in step 4) is as follows:

[0041] 4.1) The viscosity is limited by the yield stress, gradually decreasing as the yield stress increases. When the yield stress reaches...

[0042] At the limit of yield stress, the viscosity decreases as follows:

[0043]

[0044] Among them, Nd q P represents the yield viscosity. q Indicates yield stress. The second invariant of the strain rate tensor is represented; this mechanism helps to spontaneously generate weakened regions similar to real plate boundaries in the simulation, thereby enabling local plate fracture, slip, and the formation of plate boundaries.

[0045] Ultimately, in the numerical solution, the minimum value among the baseline viscosity, viscosity jump, and yield viscosity is taken as the effective viscosity of each cell.

[0046] Nd eff(d,Wj)=min(Nd(d,Wj),Nd jump (d,Wj),Nd q (9);

[0047] Therefore, no matter how large the viscosity jump is, it will not be masked by the subsequent plastic softening; at the same time, in pseudoplastic rheology, when the stress exceeds the yield threshold, the local viscosity decreases rapidly, but it does not affect the large-scale viscosity discontinuity at the phase boundary.

[0048] 4.2) Yield stress determines the critical point at which plate materials transition from elastic deformation to plastic, or even fracture behavior, and is given by the following formula:

[0049] P q =P q0 +d×dP q (10);

[0050] Among them, P q0 d represents the surface yield stress, d represents the depth, and dP q This indicates the dependence between yield stress and depth;

[0051] 4.3) Buoyancy is determined by the density difference between the tectonic plate and the surrounding mantle, and is affected by density and temperature:

[0052]

[0053] Where, m i Let m represent the density of substance i at a given surface temperature, m0 represent the reference density, and r represent the density of substance i at a given surface temperature. p ΔW represents the coefficient of thermal expansion, and ΔW represents the temperature difference in the mantle region.

[0054] The submarine plateaus with multiple parameter attributes in step 5) are divided into two categories:

[0055] 5.1) Multi-dimensional parameter attributes: Based on existing geological and geophysical exploration data, the initial length of the submarine plateau is set to 200 km, the width to 200 km, the thickness to 270 km, and the density to 3200 kg / m³, which is consistent with the density value of the subduction slab. The size of the submarine plateau is changed, defining the axis perpendicular to the subduction zone as the major axis 2x and the axis parallel to the subduction zone as the minor axis 2y, based on the initial model:

[0056] 5.1.1) Keep the minor semi-axis y unchanged, and only increase the size of the major semi-axis x;

[0057] 5.1.2) Two models with a thickness of 125 km were set up as controls;

[0058] 5.1.3) Keep the major semi-axis x unchanged, and only increase the size of the minor semi-axis y;

[0059] 5.1.4) Simultaneously increase the size of both the major and minor semi-axis;

[0060] 5.2) Multiple rheological parameter properties, focusing on three key parameters in the rheological state of seafloor plateaus: buoyancy, viscosity, and yield stress.

[0061] 5.2.1) 2900 kg / m 3 and 3200kg / m 3 As for the density variation range, with 60 kg / m³ 3 As an increment, buoyancy values ​​with different magnitudes were set according to Equation (11) to change the density of the seabed plateau;

[0062] 5.2.2) Combination of buoyancy and viscosity;

[0063] 5.2.3) Combination of buoyancy and yield stress;

[0064] 5.2.4) The combination of buoyancy, viscosity and yield stress.

[0065] The beneficial effects of this invention are:

[0066] This invention proposes and implements a method for constructing a three-dimensional spherical dynamic model of submarine plateau subduction, and provides a complete model construction and parameterization process. Through numerical comparison with traditional two-dimensional models and single-parameter three-dimensional models, this method can more comprehensively reproduce key geodynamic phenomena such as plate bending, fracturing, and arc-shaped trench formation during submarine plateau subduction. Simultaneously, it accurately captures the differentiated influences of submarine plateaus of different sizes and rheological states on local trench tectonic responses and global plate tectonics. The greatest advantage of this method lies in its ability to spontaneously generate plate tectonic behavior and plate boundary weakening through self-consistent plastic rheological approximations and multi-size, multi-rheological-state parameterized geometric-rheological scans. Furthermore, it can accurately simulate the dual impacts of submarine plateau subduction on local trenches and global plate tectonics across spatiotemporal scales. Therefore, this invention not only provides a quantitative tool for a deeper understanding of the local and global tectonic impacts of submarine plateau subduction in ocean-ocean convergence systems, but also provides solid technical support for improving plate tectonics theory, predicting regional tectonic evolution, and related seismic risk assessment and resource exploration, possessing significant theoretical value and application prospects. Attached Figure Description

[0067] Figure 1 This is a flowchart illustrating the construction principle of the three-dimensional spherical dynamic model of the submerged seabed plateau in this invention.

[0068] Figure 2 yes Figure 1 The map shows the layout of submarine plateaus of different sizes.

[0069] Figure 3It is the final constructed three-dimensional spherical dynamic model of the submerged seabed plateau.

[0070] Figure 4 This is a simulation result diagram of the impact of subduction of seafloor plateaus of different sizes on trench migration distance in an embodiment of the present invention.

[0071] Figure 5 This is a simulation result diagram of the plate fracture caused by the subduction of a large-sized seabed plateau in an embodiment of the present invention.

[0072] Figure 6 This is a comparison chart of simulation results of the impact of different buoyancy seabed plateau subduction on trench retreat in embodiments of the present invention.

[0073] Figure 7 This is a comparison chart of simulation results of the impact of subduction of seafloor plateaus with different viscosities on trench retreat in embodiments of the present invention.

[0074] Figure 8 This is a comparison chart of simulation results of the influence of different yield stresses on the subduction of the seabed plateau on the trench retreat in the embodiments of the present invention.

[0075] Figure 9 This is a simulation result diagram of the plate bending caused by the subduction of a strong rheological seafloor plateau in an embodiment of the present invention.

[0076] Figure 10 This is a simulation result diagram of the impact of the subduction of a strong rheological submarine plateau on global tectonics in an embodiment of the present invention. Detailed Implementation

[0077] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0078] A method for constructing a three-dimensional spherical dynamic model of a subducting submarine plateau includes the following steps:

[0079] Please see the flowchart of the model construction principle of this invention. Figure 1 .

[0080] 1) Solve the conservation equations for mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids;

[0081] First, the high-viscosity flow equations are solved in three-dimensional space. Specifically, the governing equations are spatially discretized on an interlaced grid, ensuring that the different unknowns to be solved are not located in the same place. Pressure and temperature are located at the center of each element volume, while the velocity components are located at the center of the corresponding element surfaces, thus avoiding artificial oscillations in the pressure field. The solver can handle relatively strong viscosity changes over small distances. To make the obtained model more realistic, viscosity can be an arbitrary function of temperature, depth, composition, and strain rate.

[0082] 2) Apply mantle convection activity

[0083] In the numerical model, the Rayleigh number is introduced as a measure of mantle convection activity.

[0084] 3) Apply viscosity jump

[0085] Viscosity can vary with temperature, depth, stress or strain rate, composition, phase transition, and melt fraction, increasing the viscosity at a depth of 660 km by 30 times.

[0086] 4) Shaping the behavior of plate-like structures

[0087] The model employs a pseudoplastic rheological approximation, enabling it to autonomously generate plate-like tectonic behavior, obtain Earth's surface velocity and tectonic activity, and form plate boundaries within and around rigid plates.

[0088] 5) Embedding submarine plateaus with multi-parameter attributes

[0089] Based on this, a global three-dimensional spherical mantle convection model was constructed. On this model, an ocean-to-ocean subduction zone was selected as the test area. A submarine plateau with multiple parameters was set up on the seaward side of the subducting ocean plate and close to the trench. This submarine plateau is characterized by multiple parameters such as geometric dimensions (length, width, and thickness) and rheological properties (buoyancy, viscosity, and yield stress).

[0090] 6) Test the model to verify the accuracy of the model solver in terms of temperature field, velocity field, buoyancy drive, etc., and finally construct a three-dimensional spherical seabed plateau subduction dynamic model;

[0091] Step 1) describes the solution of the conservation equations for mass, momentum, energy, and composition under the condition of neglecting the dynamic effects of compressible fluids, as follows:

[0092] The variable viscosity is emphasized on the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude latitude and longitude grids. The mass conservation equation, momentum conservation equation, energy conservation equation and composition conservation equation in classical fluid dynamics are solved under the condition of considering the Businesk approximation, that is, neglecting the dynamic effects of compressible fluids.

[0093] Momentum conservation equation:

[0094]

[0095] mass conservation equation:

[0096]

[0097] Energy conservation equation:

[0098]

[0099] Composition conservation equation:

[0100]

[0101] in, Let Y represent the deviatoric stress tensor, Y represent pressure, R represent the Rayleigh number (a control parameter), b represent the radius, m represent density, F represent composition, and W represent total temperature. s B represents speed. r Ha represents specific heat capacity. b Surface dissipation number, r p The coefficient of thermal expansion, s b Represents radial velocity, r c N represents the thermal conductivity coefficient, and N represents the internal heating rate. Denotes the strain rate tensor, Δm the This represents the fractional density as a function of temperature.

[0102] Step 2) describes the application of mantle convection activity as follows:

[0103] The Rayleigh number is used to characterize the convective activity of the mantle in a model, and is defined as follows:

[0104]

[0105] Where G represents gravitational acceleration, m0 represents the reference density, and r p0 The coefficient of thermal expansion of the surface is represented by ΔW, where ΔW represents the temperature difference of the mantle region of thickness M, M represents the thickness of the mantle region, and r is the temperature coefficient of thermal expansion of the surface. k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

[0106] Step 3) describes the application of viscosity jump as follows:

[0107] 3.1) Baseline viscosity: Viscosity varies radially and laterally, increases exponentially with depth, and is also temperature-dependent, defined as follows:

[0108]

[0109] Where d represents depth, Wj represents absolute temperature, and A e represents the activation energy, and Q represents the gas constant.

[0110] 3.2) Viscosity jump: At the 660km phase boundary, the baseline viscosity is multiplied by the jump factor f. jump (e.g., 30×) to reflect the viscosity discontinuity caused by rock phase transformation. The viscosity with jump correction is then obtained:

[0111] Nd jump (d,Wj)=Nd(d,Wj)fjjump (7).

[0112] The specific behavior of shaping the plate in step 4) is as follows:

[0113] 4.1) The viscosity is limited by the yield stress, gradually decreasing as the yield stress increases. When the yield stress reaches...

[0114] At the limit of yield stress, the viscosity decreases as follows:

[0115]

[0116] Among them, Nd q P represents the yield viscosity. q Indicates yield stress. This represents the second invariant of the strain rate tensor; this mechanism helps to spontaneously generate weakened regions similar to real plate boundaries in simulations, thereby enabling localized plate fracture, slip, and the formation of plate boundaries.

[0117] Ultimately, in the numerical solution, the minimum value among the baseline viscosity, viscosity jump, and yield viscosity is taken as the effective viscosity of each cell.

[0118] Nd eff (d,Wj)=min(Nd(d,Wj),Nd jump (d,Wj),Nd q (9);

[0119] Therefore, no matter how large the viscosity jump is, it will not be masked by subsequent plastic softening; at the same time, in pseudoplastic rheology, when the stress exceeds the yield threshold, the local viscosity decreases rapidly, but it does not affect the large-scale viscosity discontinuity at the phase boundary.

[0120] 4.2) Yield stress determines the critical point at which plate materials transition from elastic deformation to plastic, or even fracture behavior, and is given by the following formula:

[0121] P q =P q0 +d×dP q (10);

[0122] Among them, P q0 d represents the surface yield stress, d represents the depth, and dP q This indicates the dependence between yield stress and depth.

[0123] 4.3) Buoyancy is determined by the density difference between the tectonic plate and the surrounding mantle, and is affected by density and temperature:

[0124]

[0125] Where, m i Let m represent the density of substance i at a given surface temperature, m0 represent the reference density, and r represent the density of substance i at a given surface temperature. p ΔW represents the coefficient of thermal expansion, and ΔW represents the temperature difference in the mantle region.

[0126] The multi-parameter submarine plateau in step 5) can be specifically divided into two categories:

[0127] 5.1) Multi-dimensional parameter class: Based on existing geological and geophysical exploration data, the initial length of the submarine plateau is set to 200 km, the width to 200 km, the thickness to 270 km, and the density to 3200 kg / m³. 3 This density is consistent with the density value of the subduction slab; the first set of models, serving as a control model, does not include any seafloor plateau, so the buoyancy value of the seafloor plateau is set to 0, resulting in the control model M0. The second set of models mainly changes the size (length, width, and thickness) of the seafloor plateau, defining the axis perpendicular to the subduction zone as the major axis (2x) and the axis parallel to the subduction zone as the minor axis (2y), based on the initial model (see detailed seafloor plateau size settings for...). Figure 2 ):

[0128] 5.1.1) Keep the minor semi-axis y unchanged, and only increase the size of the major semi-axis x to obtain model M2;

[0129] 5.1.2) Two models (M1t and M2t) with a relatively thin thickness (125km) were set up as controls;

[0130] 5.1.3) Keep the major semi-axis x unchanged, and only increase the size of the minor semi-axis y to obtain model M5;

[0131] 5.1.4) Simultaneously increase the size of the major and minor semi-axis to obtain model M3;

[0132] 5.1.5) Based on model M3, keep the minor semi-axis y unchanged and only increase the size of the major semi-axis x to obtain model M4;

[0133] 5.1.6) Keep the major semi-axis x unchanged, and only increase the size of the minor semi-axis y to obtain model M6;

[0134] 5.1.7) Based on model M6, increase the size of the major and minor semi-axis to obtain model M7;

[0135] 5.2) Multiple rheological parameter class, mainly targeting three key parameters (buoyancy, viscosity, and yield stress) in the rheological state related to seafloor plateaus:

[0136] 5.2.1) 2900 kg / m 3 and 3200kg / m 3As for the density variation range, with 60 kg / m³ 3 As an increment, according to Equation (9), the density of the seabed plateau was changed, and a total of 6 models with different buoyancy values ​​were set up: M2tb-1.3, M2tb-1.1, M2tb-0.9, M2tb-0.7, M2tb-0.5, and M2t.

[0137] 5.2.2) Combination of buoyancy and viscosity (model M2tv10);

[0138] 5.2.3) Combination of buoyancy and yield stress (model M2b0s10);

[0139] 5.2.4) Combination of buoyancy, viscosity, and yield stress (model M2s10v10);

[0140] After the above steps, a three-dimensional spherical dynamic model of the subduction of the submarine plateau can be obtained (see...). Figure 3 ).

[0141] Example 1

[0142] To verify the effectiveness and correctness of "a method for constructing a three-dimensional spherical dynamic model of submarine plateau subduction," experiments and comparisons were conducted on the impacts of submarine plateau subduction on local trench migration, trench retreat, plate bending, plate fracturing, and global plate tectonics. The specific process is as follows:

[0143] 1) Subduction analysis of submarine plateaus of different sizes (see simulation results for details) Figure 4 Model M0 serves as a control model, without any submarine plateau subduction (buoyancy set to Fu = 0). Over time, the trench recedes, with the furthest distance being approximately 1050 ± 50 km (30 Myr).

[0144] 1.1) Compared with model M0, the trench retreat process in all models was greatly delayed. Different sizes of submarine plateaus subducted and corresponded to different trench retreat distances. For models M4, M6 and M7 with larger submarine plateau sizes, the trenches were delayed by about 450±50km, 850±50km and 1050±50km respectively at 30Myr.

[0145] 1.2) Models M1t to M4 exhibit similar monotonically increasing patterns. In M4, the increase in trench retreat distance gradually slows after 15 Myr, while the large-scale submarine plateau models M6 and M7 initially experience trench retreat, followed by some degree of advancement towards the trench after 22 Myr and 9 Myr, respectively. The submarine plateaus in models M1t to M3 have completely subducted before 20 Myr, while the submarine plateaus in the large-scale models M4, M6, and M7 have not completely subducted within 30 Myr.

[0146] 1.3) Model M2, with its longer seafloor plateau, has a smaller retraction distance than model M1, while model M3, with its wider seafloor plateau, has a smaller retraction distance than model M2. In models M5 and M6, the trench retraction distance is consistently smaller than that of models M2 and M4, indicating that the width of the seafloor plateau has a greater impact on its subduction effect than its length. Models M1t and M2t, with their thinner seafloor plateaus, have larger retraction distances than models M1 and M2, suggesting that thicker seafloor plateaus tend to delay trench retraction.

[0147] 1.4) Model M7, which features a large-scale submerged sea plateau, may experience plate breakage (see simulation results for details). Figure 5 ).

[0148] 2) Subduction analysis of seafloor plateaus under different rheological conditions: Buoyancy, viscosity and yield stress were incorporated into the model M2 / M2t to study the subduction of seafloor plateaus under different rheological conditions.

[0149] 2.1) First, let's consider the effect of a single buoyancy force (see simulation results for details). Figure 6 The submerged low-buoyancy seafloor plateaus (Fu = 0.0, -0.3, and -0.5) do not affect the trench geometry during their subduction process; while the high-buoyancy seafloor plateaus (Fu = -0.7, -0.9, -1.1, and -1.3) form local arc-shaped trenches in the direction of advancement, enclosing the edges of the seafloor plateaus; as the buoyancy increases, the curvature of the trench gradually increases, while the geometry of the remaining part of the trench remains unchanged.

[0150] 2.2) When the buoyancy of the seafloor plateau is small, it is insufficient to cause significant changes in the geometry of the trench (see simulation results for details). Figure 7 Using a viscosity 10 times greater to enhance the strength of the seabed plateau, simulation results show that an arc-shaped trench begins to form on the front side of the seabed plateau in the direction of subduction.

[0151] 2.3) Maintaining zero buoyancy at the seabed plateau, a dive was performed with a yield stress of 10 times the normal value (see simulation results for details). Figure 8 The results show that the trench rapidly takes on an arc-shaped form concave towards the ocean side, while the trench shape does not undergo significant deformation in the portion far from the seafloor plateau. Furthermore, the change in trench geometry is essentially a manifestation on the surface of deep-seated slab deformation; the subducting slab also exhibits relatively greater deformation and slab bending (see simulation results for details). Figure 9 ).

[0152] 3) Potential impacts on global plate tectonics (see simulation results for details) Figure 10A strong rheological model, M2s10v10, incorporating buoyancy, viscosity, and yield stress, was established and run for 200 Myr to investigate how submarine plateaus influence global plate tectonics under strong rheological conditions. After 100 Myr, subduction zone termination and mid-ocean ridge disappearance began to occur on the opposite side of the submarine plateau subduction location (opposite to the spherical view). After 200 Myr, in addition to mid-ocean ridges and subduction zone disappearance, new mid-ocean ridges and subduction zones also formed. This series of events indicates that submarine plateau subduction events occur locally and then propagate to other regions, significantly influencing plate boundary processes even far from the subduction event, and this influence requires at least 100 Myr to occur. On the other hand, this plate tectonics influence can also serve as a significant factor leading to the bending of large global seamount chains, such as the bending morphology of the Hawaii-Emperor Seamount Chain.

[0153] Therefore, the method of this invention, in constructing a self-consistent three-dimensional geodynamic model with plate-like tectonic behavior, fully considers the multi-scale geometric features and various rheological state parameters of submarine plateaus. It can accurately reproduce the tectonic response mechanism during subduction while simulating dynamic phenomena such as plate fracturing and trench bending. This method not only fills the gap in existing models' insufficient description of plateau subduction processes in ocean-ocean convergence systems but also provides a quantitative tool for a deeper understanding of the coupling mechanism between regional tectonic deformation and global plate reorganization, possessing significant theoretical value and promising prospects for widespread application.

[0154] The technical features of the above embodiments can be further combined. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

Claims

1. A method for constructing a three-dimensional spherical dynamic model of a subducting seabed plateau, characterized in that, The steps are as follows: 1) Solve the conservation equations for mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids; First, the high-viscosity flow equations are solved in three-dimensional space, that is, the governing equations are spatially discretized on an interlaced grid, and the different unknowns to be solved are not in the same location; pressure and temperature are at the center of each element volume, while velocity components are at the center of the corresponding element surface, avoiding artificial oscillations in the pressure field; the solver can handle fairly strong viscosity changes over small distances, and to make the obtained model more realistic, viscosity is an arbitrary function of temperature, depth, composition, and strain rate; 2) Apply mantle convection activity In the numerical model, the Rayleigh number is introduced as a measure of mantle convection activity; 3) Apply viscosity jump Viscosity can vary with temperature, depth, stress or strain rate, composition, phase transition, and melt fraction, increasing the viscosity at a depth of 660 km by 30 times; 4) Shaping the behavior of plate-like structures The model employs a pseudoplastic rheological approximation, enabling it to autonomously generate plate-like tectonic behavior, obtain Earth's surface velocity and tectonic activity, and form plate boundaries within and around rigid plates. The shaping behavior of the plate described in step 4) is as follows: 4.1) The viscosity is limited by the yield stress, gradually decreasing as the yield stress increases. When the yield stress reaches its limit, the viscosity decreases as follows: Among them, Nd q P represents the yield viscosity. q p represents the yield stress. ∏ The second invariant of the strain rate tensor is represented; this mechanism helps to spontaneously generate weakened regions similar to real plate boundaries in the simulation, thereby enabling local plate fracture, slip, and the formation of plate boundaries. Ultimately, in the numerical solution, the minimum value among the baseline viscosity, viscosity jump, and yield viscosity is taken as the effective viscosity of each cell. Nd eff (d,Wj)=min(Nd(d,Wj),Nd jump (d,Wj),Nd q ) (9); Therefore, no matter how large the viscosity jump is, it will not be masked by the subsequent plastic softening; at the same time, in pseudoplastic rheology, when the stress exceeds the yield threshold, the local viscosity decreases rapidly, but it does not affect the large-scale viscosity discontinuity at the phase boundary. 4.2) Yield stress determines the critical point at which plate materials transition from elastic deformation to plastic, or even fracture behavior, and is given by the following formula: P q =P q0 +d×dP q (10); Among them, P q0 d represents the surface yield stress, d represents the depth, and dP q This indicates the dependence between yield stress and depth; 4.3) Buoyancy is determined by the density difference between the tectonic plate and the surrounding mantle, and is affected by density and temperature: Where, m i Let m represent the density of substance i at a given surface temperature, m0 represent the reference density, and r represent the density of substance i at a given surface temperature. p ΔW represents the coefficient of thermal expansion, and ΔW represents the temperature difference in the mantle region. 5) Embedding submarine plateaus with multi-parameter attributes Based on this, a global three-dimensional spherical mantle convection model was constructed. On the basis of this model, an ocean-to-ocean subduction zone was selected as the test area. A submarine plateau with multiple parameters was set up on the seaward side of the subducting ocean plate and close to the trench. That is, the submarine plateau is characterized by multiple parameters including geometric dimensions such as length, width and thickness and rheological properties including buoyancy, viscosity and yield stress. 6) Test the model to verify the accuracy of the model solver in terms of temperature field, velocity field, buoyancy drive, etc., and finally construct a three-dimensional spherical seabed plateau subduction dynamic model.

2. The method according to claim 1, characterized in that, Step 1) involves solving the conservation equations for mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids. The specific steps are as follows: The variable viscosity is emphasized on the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude latitude and longitude grids. The mass conservation equation, momentum conservation equation, energy conservation equation and composition conservation equation in classical fluid dynamics are solved under the condition of considering the Businesk approximation, that is, neglecting the dynamic effects of compressible fluids. Momentum conservation equation: mass conservation equation: Energy conservation equation: Composition conservation equation: in, Let Y represent the deviatoric stress tensor, Y represent pressure, R represent the Rayleigh number (a control parameter), b represent the radius, m represent density, F represent composition, and W represent total temperature. s B represents speed. r Ha represents specific heat capacity. b Surface dissipation number, r p The coefficient of thermal expansion, s b Represents radial velocity, r c N represents the thermal conductivity coefficient, and N represents the internal heating rate. Denotes the strain rate tensor, Δm the This represents the fractional density as a function of temperature.

3. The method according to claim 1, characterized in that, Step 2) describes the application of mantle convection activity as follows: The Rayleigh number is used to characterize the convective activity of the mantle in a model, and is defined as follows: Where G represents gravitational acceleration, m0 represents the reference density, and r p0 The coefficient of thermal expansion of the surface is represented by ΔW, where ΔW represents the temperature difference of the mantle region of thickness M, M represents the thickness of the mantle region, and r represents the surface temperature. k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

4. The method according to claim 1, characterized in that, Step 3) describes the application of viscosity jump as follows: 3.1) Baseline viscosity: Viscosity varies radially and laterally, increases exponentially with depth, and is also temperature-dependent, defined as follows: Where d represents depth, Wj represents absolute temperature, and A e The activation energy is represented by Q, and the gas constant is represented by Q. 3.2) Viscosity jump: At the 660km phase boundary, the baseline viscosity is multiplied by the jump factor f. jump This reflects the viscosity discontinuity caused by the rock phase transformation; the viscosity with jump correction is then obtained: Nd jump (d,Wj)=Nd(d,Wj)f jump (7)。 5. The method according to claim 1, characterized in that, The submarine plateaus with multiple parameter attributes in step 5) are divided into two categories: 5.1) Multi-dimensional parameter attributes: Based on existing geological and geophysical exploration data, the initial length of the submarine plateau is set to 200 km, the width to 200 km, the thickness to 270 km, and the density to 3200 kg / m³, which is consistent with the density value of the subduction slab. The size of the submarine plateau is changed, defining the axis perpendicular to the subduction zone as the major axis 2x and the axis parallel to the subduction zone as the minor axis 2y, based on the initial model: 5.1.1) Keep the minor semi-axis y unchanged, and only increase the size of the major semi-axis x; 5.1.2) Two models with a thickness of 125 km were set up as controls; 5.1.3) Keep the major semi-axis x unchanged, and only increase the size of the minor semi-axis y; 5.1.4) Simultaneously increase the size of both the major and minor semi-axis; 5.2) Multiple rheological parameter properties, focusing on three key parameters in the rheological state of seafloor plateaus: buoyancy, viscosity, and yield stress. 5.2.1) 2900 kg / m 3 and 3200kg / m 3 As for the density variation range, with 60 kg / m³ 3 As an increment, buoyancy values ​​with different magnitudes were set according to Equation (11) to change the density of the seabed plateau; 5.2.2) Combination of buoyancy and viscosity; 5.2.3) Combination of buoyancy and yield stress; 5.2.4) The combination of buoyancy, viscosity and yield stress.

Images