## 3 Methods
### Ore-Forming Simulation  
#### 3.3.1 Numerical Method  
Typically, the complex ore-forming process is considered to be thermal-hydraulic-mechanical-chemical coupling processes [37; 44]. In a numerical simulation study, such processes can be approximated and characterized by a series of partial differential equations, including Darcy's Law (describing fluid activities), Fourier's Law (describing thermal convection), the Mohr-Coulomb constitutive model (describing the mechanical behavior of materials), and the mass and energy conservation equation (describing the coupled effects of the thermal-hydraulic-mechanical and chemical processes). These governing equations have been specifically reviewed in many prior studies [36; 38; 41; 44; 73; 74; 75; 76], and further elaboration is beyond the scope of the study. The basic principle of numerical simulation is to discretize the geological body model of the study area into grids or elements. At each defined grid point, the partial differential equations are transformed into algebraic equations using the finite element method or the finite difference method. These equations are subsequently solved using a number of variable parameters as initial and boundary conditions (e.g., hydraulic potential, temperature, and pore pressure; [36; 76; 77]). This process can be carried out in a series of codes or commercial software.  
Given the impracticality of universally simulating all relevant natural processes, it is necessary to choose the most appropriate software or codes related to numerical simulation, for specifically defined geological processes. In this study, the FLAC\({}^{\text{3D}}\) 5.0 software (Fast Lagrangian Analysis of Continua in Three Dimensions [78], was employed. This powerful software allows for the simulation of fluid flow and heat transport in relation to tectonic deformation, and has been successfully applied to simulation in hydrothermal-structural coupling ore-forming systems [38; 73; 79]. FLAC\({}^{\text{3D}}\) employs hexahedral elements for mesh discretization, which ensures good convergence and stability in the numerical calculations. In some scenarios with more complex geometries, polyhedral zones of other shapes (e.g., tetrahedral-shaped and pyramid zones) are deployed to create models. These discretizations can be performed in other software platforms such as SKUA-GOCAD 17 or ANSYS, and subsequently converted into a data format suitable for viewing and processing in FLAC\({}^{\text{3D}}\) by computer programs. The mechanical, hydraulic and thermodynamic parameters of the materials such as density, permeability, porosity and conductivity are assigned to the zone, representing the nature of rocks in the actual ore-forming system. The simulation process in FLAC\({}^{\text{3D}}\) proceeds through a cyclic iterative approach. During each iteration, FLAC\({}^{\text{3D}}\) computes stress, strain, and other physical quantities for each element based on the selected constitutive model, initial and boundary conditions. It then updates nodal  
\begin{table}
\begin{tabular}{l l l l} \hline \hline
**Variable** & **Definitions** & **Unit** & **Range/Mean** \\ \hline Au & Gold grade of one voxel & g/t & 0–13.4/0.92 \\ AuMet & Gold amounts of one voxel & g & 0–303,532/7784 \\ \(d\)F & Minimum distance to the Zhaoping fault & m & \(-416\)–113/ \(-117\) \\ \(g\)F & Dip of the Zhaoping fault & \({}^{\circ}\) & 8.8–89/47 \\ \(v\)F & Dip variation of the Zhaoping fault in a 100 m buffer & \({}^{\circ}\)/100 m & \(-15.2\)–45.6/3.2 \\ \(wa\)F & Undulation relative to trending surface of the Zhaoping fault in a 180 m buffer & m & \(-164\)–96/ \(-15.7\) \\ \(wb\)F & Undulation relative to trending surface of the Zhaoping fault in a 360 m buffer & m & \(-111\)–\(101/-13.5\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Model geometric feature and ore attribute description.  
displacements and other variables according to results obtained from the execution of the previous cycle [80].