Method for constructing three-dimensional spherical dynamic model of seabed plateau dive

Abstract

The invention discloses a method for constructing a three-dimensional spherical dynamic model of seabed plateau dive. The method comprises the following steps: solving mass, momentum, energy and component conservation equations under the condition of neglecting the dynamic influence of compressible fluid; applying mantle convection vitality; applying a viscosity jump; an autonomous plate-like construction behavior is shaped; embedding a seabed plateau with multi-parameter attributes; and finally completing the construction of the three-dimensional spherical earth dynamics model after model testing. The constructed three-dimensional numerical model adopts a three-dimensional spherical construction method with an autonomous class plate construction behavior, and fine description is performed on the seabed plateau by introducing different sizes and different rheological state parameters. The model not only can accurately reproduce the geodynamic phenomena such as bending and fracture of the plate in the subduction process of the seabed plateau, but also can reveal the response difference of the seabed plateau with different sizes and rheological characteristics to the local sea ditch structure, and realizes the cross-scale simulation of the influence from the local structure to the global structure.

Description

Technical Field

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

[0002] Given the difficulty and high cost of directly acquiring deep Earth information, numerical simulation has long been the mainstream international approach for understanding the dynamics of submarine plateau subduction. This type of simulation not only reproduces the dynamic evolution and deformation characteristics revealed by observational and survey data, but also transcends the limitations of geological time and space, enabling comprehensive exploration of a wide range of Earth science issues and providing quantitative support for observational results, making it an indispensable tool in marine geoscience research. Numerical simulations are generally categorized into two types based on spatial dimension: two-dimensional and three-dimensional.

[0003] Numerous studies have focused on the subduction of oceanic plates carrying submarine plateaus using two-dimensional models. Yan et al. (2021) used two-dimensional thermodynamic simulations to find that when a submarine plateau possesses a thick lithospheric basement and the continental plate undergoes thrusting, the insufficient buoyancy of the plateau can hinder subduction and ultimately trigger a jump. Dai et al. (2020) also used a two-dimensional thermodynamic model to demonstrate that large submarine plateaus on the Paleo-Pacific Plate could induce flat subduction, which evolved into steep subduction from the Mesozoic to the Cenozoic. However, because such simulations assume no change in trench direction, they highly simplify the geological process and thus have limited applicability. In contrast, three-dimensional numerical simulations, with their greater complexity and sophistication, can produce results that are more realistic and authentic to nature. For example, Liu et al. (2010), combining three-dimensional simulations with paleoplate reconstructions, found that the timing of the subduction and removal of the Shatsky submarine plateau coincided with the Laramide orogeny, suggesting that plateau subduction could have triggered regional surface rebound and uplift. Betts et al. (2015) introduced the mantle plume effect in three-dimensional experiments, further revealing the key influence of the synergistic effect of submarine plateaus and mantle plumes on the geometric evolution of trenches and subducting plates.

[0004] However, current two-dimensional and three-dimensional studies of the effects of submarine plateau subduction have largely focused on ocean-continent convergence zones, with insufficient attention paid to processes underlying ocean-ocean convergence. Furthermore, existing three-dimensional models often neglect the role 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-scale, multi-rheological submarine plateaus and subduction zones has not been systematically explored. Consequently, existing three-dimensional simulations still fail to fully capture the overall impact of submarine plateau subduction on the tectonic and dynamic behavior of subduction zones. Summary of the Invention

[0005] This paper discloses a method for constructing a three-dimensional spherical dynamic model of submarine plateau subduction, aiming to address the problem that existing three-dimensional model studies ignore the convergence of oceanic plates and the effects of overlying plates. By using viscosity jumps and pseudoplastic rheological approximations, a self-consistent plate-like tectonic behavior is generated. Furthermore, submarine plateaus with multiple sizes and rheological parameters are introduced to achieve the construction of a three-dimensional spherical dynamic model of submarine plateau subduction.

[0006] The present invention is achieved through the following technical solutions:

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

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

[0009] First, the high-viscosity flow equations are solved in three dimensions. Specifically, the governing equations are spatially discretized on a staggered grid. The different unknowns to be solved are not located in the same location. The pressure and temperature are located at the center of each unit volume, while the velocity components are located at the center of the corresponding unit face. This avoids the appearance of artificial oscillations in the pressure field. The solver can handle quite strong viscosity variations over small distances. To make the resulting model more realistic, the viscosity can be an arbitrary function of temperature, depth, composition, and strain rate.

[0010] 2) Imposing mantle convection activity

[0011] In numerical models, 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 change, and melt fraction, increasing the viscosity at a depth of 660 km by a factor of 30;

[0014] 4) Shaping plate-like structural behavior

[0015] The model uses a pseudoplastic rheological approximation, which enables the model to autonomously generate plate-like tectonic behavior, obtain the Earth's surface velocity and tectonic action, 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 the basis of this model, an ocean-to-ocean subduction zone was selected as the experimental area. A submarine plateau with multi-parameter properties was set up on the seaward side of the subducting oceanic plate and close to the trench. That is, the submarine plateau was comprehensively characterized by multiple parameters such as geometric dimensions (length, width, thickness) and rheological properties (buoyancy, viscosity, 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 submarine plateau subduction dynamics model.

[0019] The mass, momentum, energy and composition conservation equations described in step 1) are solved while ignoring the dynamic effects of the compressible fluid as follows:

[0020] Emphasis is placed on the variable viscosity of the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude longitude and latitude grids. The Boussinesq approximation, which ignores the dynamic effects of compressible fluids, is used to solve the mass conservation equations, momentum conservation equations, energy conservation equations, and composition conservation equations in classical fluid dynamics.

[0021] Momentum conservation equation:

[0022]

[0023] The mass conservation equation:

[0024]

[0025] Energy conservation equation:

[0026]

[0027] Composition conservation equation:

[0028]

[0029] in, represents the deviatoric stress tensor, Y represents pressure, R represents the controlling parameter Rayleigh number, b represents radius, m represents density, F represents composition, W represents total temperature, s Indicates speed, B r Represents specific heat capacity, Ha b Surface dissipation number, r p represents the thermal expansion coefficient, s b represents the radial velocity, r c represents the thermal conductivity, N represents the internal heating rate, represents the strain rate tensor, Δm the represents the fractional density that varies with temperature,

[0030] The step 2) of applying mantle convection energy is specifically as follows:

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

[0032]

[0033] Where G represents the acceleration of gravity, m0 represents the reference density, r p0 represents the thermal expansion coefficient of the surface, ΔW represents the temperature difference of the mantle region with a thickness of M, M represents the thickness of the mantle region, r k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

[0034] The viscosity jump applied in step 3) is as follows:

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

[0036]

[0037] Where d represents the depth, Wj represents the absolute temperature, A e represents activation energy, Q represents gas constant;

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

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

[0040] The plate-like structural behavior described in step 4) is as follows:

[0041] 4.1) The viscosity is limited by the yield stress and gradually decreases with the increase of the yield stress. When the yield stress reaches

[0042] When the yield stress reaches its limit, the viscosity decreases to:

[0043]

[0044] Among them, Nd q represents the yield viscosity, P q represents the yield stress, represents the second invariant of the strain rate tensor; this mechanism helps to spontaneously generate weakened regions similar to real plate boundaries in the simulation, thereby achieving localized plate breakage, slip, and the formation of plate boundaries;

[0045] Finally, 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, when the stress exceeds the yield threshold in pseudoplastic rheology, the local viscosity decreases rapidly without affecting the large-scale viscosity discontinuity at the phase boundary;

[0048] 4.2) The yield stress determines the critical point at which the plate material transitions from elastic deformation to plasticity, or even fracture behavior, and is given by the following equation:

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

[0050] Among them, P q0 represents the yield stress of the surface, d represents the depth, dP q represents the dependence of yield stress on depth;

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

[0052]

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

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

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

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

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

[0058] 5.1.3) Keep the major axis x unchanged and only increase the minor axis y;

[0059] 5.1.4) Increase the size of both the major and minor semi-axes simultaneously;

[0060] 5.2) Multi-rheological parameter attributes, targeting the three key parameters of buoyancy, viscosity, and yield stress in the rheological state related to submarine plateaus:

[0061] 5.2.1) 2900kg / m 3 and 3200kg / m 3 As the density variation range, 60kg / m 3 For the increment, different buoyancy values are set according to formula (11) by changing the density of the submarine plateau;

[0062] 5.2.2) Combination of buoyancy and viscosity;

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

[0064] 5.2.4) Combination of buoyancy with viscosity and yield stress.

[0065] Beneficial effects of the present invention:

[0066] This paper 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, fracture, and arc trench formation during submarine plateau subduction. It also accurately captures the differentiated impacts of submarine plateaus of different sizes and rheological states on local trench tectonic responses and global plate tectonics. The greatest advantage of this method is that, through a self-consistent plastic rheological approximation and multi-size, multi-rheological state parameterized geometric-rheological scanning, it can not only spontaneously generate plate-like tectonic behavior and plate boundary weakening, but also accurately simulate the dual effects of submarine plateau subduction on local trenches and global plate tectonics across time and space scales. Therefore, this method not only provides a quantitative tool for in-depth 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 tectonic theory, predicting regional tectonic evolution, related earthquake risk assessment, and resource exploration. It has important theoretical value and application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 This is a flow chart of the principle of constructing a three-dimensional spherical dynamic model of submarine plateau subduction according to the present invention.

[0068] Figure 2 yes Figure 1 A diagram showing the settings of submarine plateaus of different sizes.

[0069] Figure 3This is the final constructed three-dimensional spherical dynamic model of submarine plateau subduction.

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

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

[0072] Figure 6 This is a comparison chart of simulation results of the influence of subduction of submarine plateaus with different buoyancy on trench retreat in an embodiment of the present invention.

[0073] Figure 7 This is a comparison chart of simulation results of the influence of submarine plateau subduction on trench retreat at different viscosities in an embodiment of the present invention.

[0074] Figure 8 This is a comparison chart of simulation results of trench retreat affected by subduction of submarine plateaus with different yield stresses in an embodiment of the present invention.

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

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

[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 submarine plateau subduction includes the following steps:

[0079] The flow chart of the model construction principle of the present invention is shown in Figure 1 .

[0080] 1) Solve the conservation equations of 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 a staggered grid. The different unknowns to be solved are not located in the same location. The pressure and temperature are located at the center of each unit volume, while the velocity components are located at the center of the corresponding unit face. This avoids artificial oscillations in the pressure field. The solver can handle quite strong viscosity changes over small distances. To make the resulting model closer to reality, the viscosity can be an arbitrary function of temperature, depth, composition, and strain rate.

[0082] 2) Imposing mantle convection activity

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

[0084] 3) Apply viscosity jump

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

[0086] 4) Shaping plate-like structural behavior

[0087] The model uses a pseudoplastic rheological approximation, which enables the model to autonomously generate plate-like tectonic behavior, obtain the Earth's surface velocity and tectonic action, 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 the basis of this model, an ocean-to-ocean subduction zone was selected as the experimental area. A submarine plateau with multi-parameter properties was set up on the seaward side of the subducting oceanic plate and close to the trench. That is, the submarine plateau was comprehensively characterized by multiple parameters such as geometric dimensions (length, width, thickness) and rheological properties (buoyancy, viscosity, 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 submarine plateau subduction dynamics model;

[0091] The solution of the mass, momentum, energy and composition conservation equations described in step 1) while ignoring the dynamic effects of the compressible fluid is as follows:

[0092] Emphasis is placed on the variable viscosity of the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude longitude and latitude grids. The Boussinesq approximation, which ignores the dynamic effects of compressible fluids, is used to solve the mass conservation equations, momentum conservation equations, energy conservation equations, and composition conservation equations in classical fluid dynamics.

[0093] Momentum conservation equation:

[0094]

[0095] The mass conservation equation:

[0096]

[0097] Energy conservation equation:

[0098]

[0099] Composition conservation equation:

[0100]

[0101] in, represents the deviatoric stress tensor, Y represents pressure, R represents the controlling parameter Rayleigh number, b represents radius, m represents density, F represents composition, W represents total temperature, s Indicates speed, B r Represents specific heat capacity, Ha b Surface dissipation number, r p represents the thermal expansion coefficient, s b represents the radial velocity, r c represents the thermal conductivity, N represents the internal heating rate, represents the strain rate tensor, Δm the represents the fractional density that varies with temperature,

[0102] The step 2) of applying mantle convection energy is specifically as follows:

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

[0104]

[0105] Where G represents the acceleration of gravity, m0 represents the reference density, r p0 represents the thermal expansion coefficient of the surface, ΔW represents the temperature difference of the mantle region with a thickness of M, M represents the thickness of the mantle region, r k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

[0106] The viscosity jump described in step 3) is specifically as follows:

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

[0108]

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

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

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

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

[0113] 4.1) The viscosity is limited by the yield stress and gradually decreases with the increase of the yield stress. When the yield stress reaches

[0114] When the yield stress reaches its limit, the viscosity decreases to:

[0115]

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

[0117] Finally, 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 the subsequent plastic softening; at the same time, when the stress exceeds the yield threshold in pseudoplastic rheology, the local viscosity decreases rapidly without affecting the large-scale viscosity discontinuity at the phase boundary.

[0120] 4.2) The yield stress determines the critical point at which the plate material transitions from elastic deformation to plasticity, or even fracture behavior, and is given by the following equation:

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

[0122] Among them, P q0 represents the yield stress of the surface, d represents the depth, dP q Represents the dependence of yield stress on depth.

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

[0124]

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

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

[0127] 5.1) Multi-scale parameter class: Based on existing geological and geophysical exploration data, the initial length of the submarine plateau is set to 200 km, the width is set to 200 km, the thickness is set to 270 km, and the density is set to 3200 kg / m 3 , which is consistent with the density of the subducting plate. The first set of models is used as a control model. No submarine plateau is placed, so the buoyancy value of the submarine plateau is set to 0, and the control model M0 is obtained. The second set of models mainly changes the size (length, width and thickness) of the submarine plateau. The axis perpendicular to the subduction zone is defined as the major axis (2x), and the axis parallel to the subduction zone is defined as the minor axis (2y). Based on the initial model (for detailed submarine plateau size settings, see Figure 2 ):

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

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

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

[0131] 5.1.4) Increase the size of both the major and minor semi-axes simultaneously to obtain model M3;

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

[0133] 5.1.6) Keep the major axis x constant and only increase the minor axis y to obtain model M6;

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

[0135] 5.2) Multi-rheological parameter class, mainly focusing on three key parameters (buoyancy, viscosity and yield stress) in the rheological state related to submarine plateaus:

[0136] 5.2.1) 2900kg / m 3 and 3200kg / m 3As the density variation range, 60kg / m 3 For increment, according to formula (9), the density of submarine plateau is changed to set up 6 models with different buoyancy values, including 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, we can finally obtain a three-dimensional spherical dynamic model of submarine plateau subduction (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 effects of submarine plateau subduction on local trench migration, trench retreat, slab bending, slab fracture, and global plate tectonics. The specific process is as follows:

[0143] 1) Analysis of subduction of submarine plateaus of different sizes (see simulation results for details) Figure 4 ): Model M0 is used as the control model. There is no submarine plateau subduction (buoyancy is set to Fu=0). As time goes by, the trench retreats, and the farthest distance is about 1050±50km (30Myr).

[0144] 1.1) Compared with model M0, the trench retreat process in all models is greatly delayed. The subduction of submarine plateaus of different sizes corresponds to different trench retreat distances. For models M4, M6 and M7 with larger submarine plateaus, the trench retreat distances are delayed by approximately 450±50km, 850±50km and 1050±50km at 30Myr, respectively.

[0145] 1.2) Models M1t through M4 exhibit a similar monotonic increase. The increase in trench retreat in M4 slows after 15 Myr, while the large-scale submarine plateau models M6 and M7 initially experience trench retreat followed by a certain degree of trench advancement after 22 Myr and 9 Myr, respectively. The submarine plateaus in models M1t through M3 were fully subducted before 20 Myr, while those in large-scale models M4, M6, and M7 did not fully subduct within 30 Myr.

[0146] 1.3) Model M2, with its longer submarine plateau, has a smaller retreat distance than model M1, while model M3, with its wider submarine plateau, has a smaller retreat distance than model M2. Models M5 and M6 consistently show smaller trench retreat distances than models M2 and M4, suggesting that submarine plateau width, rather than length, has a greater influence on subduction. Models M1t and M2t, with their thinner submarine plateaus, have larger retreat distances than models M1 and M2, indicating that thicker submarine plateaus tend to delay trench retreat.

[0147] 1.4) Model M7 with a large submarine plateau subducting, which will cause the plate to break (see simulation results for details) Figure 5 ).

[0148] 2) Analysis of subduction of submarine plateaus under different rheological states: Buoyancy, viscosity and yield stress are incorporated into the model M2 / M2t to study the subduction of submarine plateaus under different rheological states.

[0149] 2.1) First, the effect of single buoyancy (see simulation results for details) Figure 6 ), weak submarine plateaus with small buoyancy (Fu=0.0, -0.3 and -0.5), their subduction process has no effect on the geometric shape of the trench; while strong submarine plateaus with higher buoyancy (Fu=-0.7, -0.9, -1.1 and -1.3) form a local arc-shaped trench shape in the forward direction, wrapping the edge of the submarine plateau; as the buoyancy continues to increase, the curvature of the trench is gradually magnified, while the geometric shape of the remaining part of the trench remains unchanged.

[0150] 2.2) When the buoyancy of the submarine plateau is small, it is not enough to cause significant changes in the geometry of the trench (see simulation results for details). Figure 7 ), using 10 times the viscosity to enhance the strength of the submarine plateau, the simulation results show that an arc-shaped trench begins to form in front of the subduction direction of the submarine plateau.

[0151] 2.3) Keep the buoyancy of the submarine plateau at 0 and set 10 times the yield stress for subduction (see simulation results for details) Figure 8 ), the results show that the trench quickly takes on an arc shape concave toward the ocean side, while the shape of the trench away from the submarine plateau does not undergo significant deformation. In addition, the change in the geometric shape of the trench is essentially due to the surface manifestation of the deformation of the deep plate. At this time, the deformation of the subducting plate is also relatively large, and the plate bending occurs (see simulation results for details). Figure 9 ).

[0152] 3) Potential impact on global plate tectonics (see simulation results for details) Figure 10): A strong rheological model M2s10v10 with buoyancy, viscosity, and yield stress was set up and run for 200 Myr to explore how submarine plateaus under strong rheological conditions affect global plate tectonics. After 100 Myr, the termination of subduction zones and the disappearance of mid-ocean ridges began to appear on the other side of the submarine plateau subduction location (the opposite side from a spherical perspective). After 200 Myr, in addition to the disappearance of ridges and subduction, new mid-ocean ridges and subduction zones were also formed. This series of events shows that submarine plateau subduction events occur in local areas and then spread to other areas. Even far away from the subduction event, they can have a significant impact on plate boundary processes, and this impact takes at least 100 Myr to occur. On the other hand, this plate tectonic influence can be an important factor in the curvature of large seamount chains around the world, such as the curved morphology of the Hawaiian-Emperor seamount chain.

[0153] This demonstrates that the proposed method fully accounts for the multiscale geometric characteristics and rheological parameters of submarine plateaus in constructing a self-consistent three-dimensional geodynamic model with plate-like tectonic behavior. This method accurately reproduces the tectonic response mechanisms during subduction while simulating dynamic phenomena such as slab fracture and trench bending. This method not only addresses the inadequate description of plateau subduction processes in ocean-ocean convergence systems by existing models but also provides a quantitative tool for a deeper understanding of the coupling mechanisms between regional tectonic deformation and global plate reorganization, possessing significant theoretical value and promising prospects for widespread application.

[0154] The various technical features of the above-described embodiments can be further combined. To make the description concise, not all possible combinations of the various technical features in the above-described 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 submarine plateau subduction, characterized in that: Here are the steps: 1) Solve the conservation equations of mass, momentum, energy, and composition while neglecting the dynamic effects of compressible fluids; First, the high-viscosity flow equations are solved in three dimensions. This involves spatially discretizing the governing equations on a staggered grid, so that the different unknowns to be solved are not located in the same location. Pressure and temperature are located at the center of each unit volume, while velocity components are located at the center of the corresponding unit face, to avoid artificial oscillations in the pressure field. The solver can handle fairly strong viscosity variations over small distances. To achieve a more realistic model, viscosity is an arbitrary function of temperature, depth, composition, and strain rate. 2) Imposing mantle convection activity In numerical models, 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 change, and melt fraction, increasing the viscosity at a depth of 660 km by a factor of 30; 4) Shaping plate-like structural behavior The model uses a pseudoplastic rheological approximation, which enables the model to autonomously generate plate-like tectonic behavior, obtain the Earth's surface velocity and tectonic action, and form plate boundaries within and around rigid plates; 5) Embedding submarine plateaus with multi-parameter attributes Based on this, a global three-dimensional spherical mantle convection model was constructed. Based on this model, an ocean-to-ocean subduction zone was selected as the experimental area. A submarine plateau with multi-parameter attributes was set up on the subducting oceanic plate and close to the seaward side of the trench. The submarine plateau was characterized by multiple parameters including geometric dimensions such as length, width, and thickness, and rheological properties such as buoyancy, viscosity, and yield stress. 6) Test the model to verify the accuracy of the model solver in terms of temperature field, velocity field, and buoyancy drive, and finally construct a three-dimensional spherical submarine plateau subduction dynamics model.

2. The method according to claim 1, characterized in that The solution of the mass, momentum, energy and composition conservation equations described in step 1) while ignoring the dynamic effects of the compressible fluid is as follows: Emphasis is placed on the variable viscosity of the yin-yang grid, which divides the sphere into two orthogonal and slightly overlapping low-latitude longitude and latitude grids. The Boussinesq approximation, which ignores the dynamic effects of compressible fluids, is used to solve the mass conservation equations, momentum conservation equations, energy conservation equations, and composition conservation equations in classical fluid dynamics. Momentum conservation equation: The mass conservation equation: Energy conservation equation: Composition conservation equation: in, represents the deviatoric stress tensor, Y represents pressure, R represents the controlling parameter Rayleigh number, b represents radius, m represents density, F represents composition, W represents total temperature, s Indicates speed, B r Represents specific heat capacity, Ha b Surface dissipation number, r p represents the thermal expansion coefficient, s b represents the radial velocity, r c represents the thermal conductivity, N represents the internal heating rate, represents the strain rate tensor, Δm the represents the fractional density that varies with temperature, 3. The method according to claim 1, characterized in that The step 2) of applying mantle convection energy is as follows: The Rayleigh number is used to characterize the convective activity of the mantle in the model and is defined as follows: Where G represents the acceleration of gravity, m0 represents the reference density, r p0 represents the thermal expansion coefficient of the surface, ΔW represents the temperature difference of the mantle region with a thickness of M, M represents the thickness of the mantle region, r k represents the thermal diffusivity, and Nd0 represents the reference viscosity.

4. The method according to claim 1, wherein The viscosity jump applied in step 3) is as follows: 3.1) Baseline Viscosity: Viscosity varies radially and laterally, increases exponentially with depth, and is temperature-dependent. It is defined as follows: Where d represents the depth, Wj represents the absolute temperature, A e represents activation energy, Q represents gas constant; 3.2) Viscosity jump: At the 660 km phase boundary, multiply the baseline viscosity by the jump factor f jump , to reflect the viscosity discontinuity caused by rock phase transition; at this time, the jump-corrected viscosity is obtained: Nd jump (d,Wj)=Nd(d,Wj)f jump (7)。 5. The method according to claim 1, wherein The plate-like structural behavior described in step 4) is as follows: 4.1) The viscosity is limited by the yield stress and gradually decreases with the increase of the yield stress. When the yield stress reaches the limit of the yield stress, the viscosity decreases to: Among them, Nd q represents the yield viscosity, P q represents the yield stress, represents the second invariant of the strain rate tensor; this mechanism helps to spontaneously generate weakened regions similar to real plate boundaries in the simulation, thereby achieving localized plate breakage, slip, and the formation of plate boundaries; Finally, 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 uump (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, when the stress exceeds the yield threshold in pseudoplastic rheology, the local viscosity decreases rapidly without affecting the large-scale viscosity discontinuity at the phase boundary; 4.2) The yield stress determines the critical point at which the plate material transitions from elastic deformation to plasticity, or even fracture behavior, and is given by the following equation: P q =P q0 +d×dP q (10); Among them, P q0 represents the yield stress of the surface, d represents the depth, dP q represents the dependence of yield stress on depth; 4.3) Buoyancy is determined by the density difference between the plate and the surrounding mantle, and is affected by density and temperature: Among them, m i represents the density of substance i at a given surface temperature, m0 represents the reference density, and r p represents the thermal expansion coefficient, and ΔW represents the temperature difference in the mantle region.

6. The method according to claim 1, characterized in that The submarine plateaus with multi-parameter attributes in step 5) are divided into two categories: 5.1) Multi-scale parameter attributes. Based on existing geological and geophysical exploration data, the initial length of the submarine plateau is set to 200 km, the width is set to 200 km, the thickness is set to 270 km, and the density is set to 3200 kg / m3, which is consistent with the density of the subducting plate. The size of the submarine plateau is changed, and the axis perpendicular to the subduction zone is defined as the major axis 2x, and the axis parallel to the subduction zone is defined as the minor axis 2y. Based on the initial model: 5.1.1) Keep the minor axis y unchanged and only increase the size of the major axis x; 5.1.2) Two models with a thickness of 125 km were set as controls; 5.1.3) Keep the major axis x unchanged and only increase the minor axis y; 5.1.4) Increase the size of both the major and minor semi-axes simultaneously; 5.2) Multi-rheological parameter attributes, targeting the three key parameters of buoyancy, viscosity, and yield stress in the rheological state related to submarine plateaus: 5.2.1) 2900kg / m 3 and 3200kg / m 3 As the density variation range, 60kg / m 3 For the increment, different buoyancy values are set according to formula (11) by changing the density of the submarine plateau; 5.2.2) Combination of buoyancy and viscosity; 5.2.3) Combination of buoyancy and yield stress; 5.2.4) Combination of buoyancy with viscosity and yield stress.

Images