To investigate the present-day stress state, we developed a numerical model based on conservation equations using COMSOL Multiphysics 5.453,54. Assuming that accelerations other than gravity can be neglected20, we solve the force-balance Eq. 55 for absolute stress (no pore pressure):

$$\begin{array}{c}\nabla \bullet \sigma +{F}_{v}=0 ,\left(1\right)\end{array}$$

where \(\sigma\) is stress tensor, and \({F}_{v}\) is body force. The compressive stress here is reckoned to be positive and the tensile stress is negative. We adopt an explicit hybrid finite element approach to solve force-balance equation. According to Hook’s Law, the relationship between the stress tensor and strain, \({\epsilon }_{ij}\), can be described as following:

$$\begin{array}{c}{\sigma }_{ij}=E{\epsilon }_{ij} ,\left(2\right)\end{array}$$

Where \(i\) and \(j\) are coordinate indices (*x, y, z*), and \(E\) is the Young’s modulus. The shear strain components can be defined as follows:

$$\begin{array}{c}{\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}\right) , \left(3\right)\end{array}$$

Where \(u\) is the displacement.

The initial numerical models comprise a 52-km-long, 53-km-wide and ~ 4.5-km-deep domain and a 13-km-long, 13-km-wide and ~ 4.5-km-deep domain, respectively (see in the supplementary materials). The geometry of the top surface in our model is consistent with the seafloor topography at a resolution of ~ 50 m. For simplification, the other surfaces are simulated as planes. The horizontal x axis is parallel to the ridge axis, the horizontal y axis is perpendicular to the ridge axis, and the vertical z axis is directed downward. In order to eliminate boundary effects, the area of interest (seafloor in our model) is set at a significant distance away from the bottom surface. The minimum and maximum thickness of our models are 2.5 and ~ 4.5 km, respectively.

Our model contains two different stratigraphic layers (Layer A and Layer B). Based on the physical properties model developed by Hyndman and Drury (1976) 56 and active-source seismic experiments23. The density of the two layers is measured from the samples 57. Yang’s modules and Poisson’s ratio are derived from active-source P-wave velocity data based on empirical equations58. The half-spreading rates are loaded on the model edges. Zero-displacement boundary condition is imposed along the x and z direction while along the y direction, both edges are slip-free21. The normal stress equals the overburden pressure and shear stress is set to be zero on the bottom surface. The detail parameters of our model are shown in Table 1 in supplementary materials. Hydrostatic pressure (P) is applied to the top surface, as follows:

$$\begin{array}{c}{P}_{i,j}={\rho }_{w}g{z}_{i,j} ,\left(4\right)\end{array}$$

where P is hydrostatic pressure, \({\rho }_{w}\) is the water density, and z is the depth. Acceleration due to gravity (\(g\)) is 9.8 \(m\bullet {\text{s}}^{-2}\).

## Data Availability

The focal mechanism data used in this study are supplied in the Supplementary Tables. The bathymetry data of TAG available on GEOMAP app. Bathymetric data of Longqi can be provided upon request by sending an e-mail to the corresponding author ([email protected]).