## 2 Materials and Methods
### Numerical Simulation
#### 2.4.2 Initial and Boundary Conditions  
This numerical simulation study focuses on the transport and transformation processes of selected reactive species (HF, H\({}_{2}\)SO\({}_{4}\)) in phosphogypsum leachate under dynamic dissolution conditions. The anisotropic flow behavior in fractured rock masses is primarily controlled by fracture network characteristics. Considering the substantial variations in fracture orientation, density, and connectivity among different rock specimens, we systematically selected specimens with nearly identical fracture attributes (including quantity, spatial distribution, and developmental direction) for dynamic dissolution experiments to ensure analytical consistency. This selection protocol justified the exclusive use of the HY10 specimen for subsequent numerical modeling. While pH variations influence reaction kinetics, they do not fundamentally modify solute transport mechanisms within the rock matrix. Consequently, the numerical simulations specifically addressed phosphogypsum leachate-carbonate rock dissolution interactions under pH = 2 conditions. The HY10 rock specimen was imaged using a multi-scale high-resolution X-ray 3D digital core imaging system (X-CT, NanoVoxel4561, Sanying Precision Instruments Co., Ltd., Tianjin, China). The reconstructed 3D fracture model generated through Avizo's image processing pipeline was subsequently exported in STL formats. These files were then imported into COMSOL Multiphysics (version 6.2) to establish a refined fracture mesh model. But the mesh model contained excessive mesh elements with incomplete curves and unclosed geometries, preventing direct conversion into a physical model. To address this, the fracture model was simplified. The rock specimen's fracture parameters, including average fracture aperture (0.50 mm), dominant fracture orientation, and cumulative fracture length, were quantified using the Label Analysis module in Avizo software (version 2020.1). The original complex fracture network was replaced with five rectangular thin plates whose apertures matched the measured average value (0.50 mm), with identical spatial orientations and positions. A cylindrical host rock measuring 85 mm in diameter and 30 mm in height, consistent with the actual specimen dimensions, was then incorporated. By varying the plate thickness while maintaining aperture consistency, three models with fracture apertures of 0.20 mm, 0.50 mm, and 0.05 mm were established. Gravitational effects were excluded from all simulations. The cylindrical rock specimen model's top surface was assigned as the inflow boundary (pressure inlet), and the bottom surface was designated as the outflow boundary (pressure outlet). Fine hierarchy in COMSOL Multiphysics was selected based on a balanced consideration of average mesh quality and solver performance, prioritizing optimal trade-offs between numerical accuracy and computational efficiency. The meshing strategy employs a physics-controlled mesh generation approach, where the software automatically generates the computational grid based on predefined boundary conditions of the Darcy flow field and solute concentration field. The model parameters were primarily derived from geotechnical experimental data, engineering reports, and empirical correlations. In COMSOL Multiphysics, the pH value and temperature were defined as global parameters, set to 2 and 293.15 K, respectively, under the Global Definitions node to establish consistent environmental conditions throughout the model. For the Darcy flow field configuration, the permeability (\(k\)), porosity (\(\epsilon_{p}\)), and specific storativity (\(\mu S\)) were adopted from Table 6, while other hydrodynamic parameters utilized standard water properties from the COMSOL Multiphysics material library, neglecting the physicochemical discrepancies between phosphogypsum leachate and deionized water. The reaction rates for chemical species must be configured within the Dilute Species Transport interface.