## 5 Discussion
One potential explanation for this discrepancy is the choice of input data representation: only eight key points from each stress-strain curve (sixteen features total) were extracted rather than utilizing full-curve information or more representative parameters. This limited sampling may omit critical nonlinear features inherent in the continuous curve, restricting the model's capacity to learn complex material behavior. In addition, all the simulation curves end at 50% of the peak. This truncation inevitably limits the model's ability to characterize \(E_{d}\), which governs behavior beyond the peak. Future work will extend simulations deeper into the post-peak regime to better capture and invert for \(E_{d}\). Additionally, we plan to pre-train models using the DEM simulation results of geomechanically similar materials (e.g., sandstone or limestone) and utilize these as warm-start initializations in the AutoML framework. This strategy aims to improve the generalization capability while reducing the training burden for new materials. Furthermore, incorporating physically meaningful derived features, such as the tangent modulus and peak slope, will be explored in future work to improve model interpretability and accuracy. Three mesoscopic parameters (\(E\), \(k_{s}/k_{n}\), and \(E_{d}\)) were selected as the model outputs in this study, whereas the linear parallel bond model in PFC3D comprises seven parameters, including the friction coefficient (\(\mu\)), tensile strength (\(\sigma t\)), cohesion (\(c\)), and post-failure stiffness ratios (\(k_{s,d}/k_{n,d}\)). Previous studies have shown that parameters such as \(\mu\), \(\sigma t\), and \(c\) also influence the macroscopic mechanical responses [29]. Omitting these parameters may therefore result in incomplete characterization of the input-output mapping, leading to suboptimal inversion performance. Moreover, nonlinear and coupled interactions between parameters, for instance, the friction coefficient and elastic modulus jointly affecting the stiffness, cannot be fully captured when only a subset is predicted. Future research should extend the inversion framework to encompass all seven mesoscopic parameters. A further consideration pertains to the black-box nature of machine learning models, including both ANNs and AutoML. Such models generally lack direct physical interpretability and may produce unrealistic predictions when deployed beyond their training domains. To address this issue, Physics-Informed Neural Networks (PINNs) offer a promis  
Figure 9: The inversion performance of the developed multi-objective ANN model on training and testing datasets: (**a**) \(E\); (**b**) \(k_{s}/k_{n}\); and (**c**) \(E_{d}\). ing direction. PINNs embed governing equations (e.g., equilibrium conditions, constitutive relationships) directly into the loss function, ensuring that network predictions satisfy known physical laws while learning from data [30]. In geomechanical contexts, PINNs have shown success in fluid mechanics, fracture modeling, and subsurface flow [31,32]. Incorporating PFC3D's constitutive framework or energy balance constraints into a PINN architecture could improve both the inversion accuracy and model robustness. By enforcing physical consistency, PINNs can also facilitate extrapolation beyond the training dataset--a capability that proves especially valuable when experimental data are limited or when simulating conditions not encompassed by the original orthogonal design [33, 34, 35, 36, 37]. Inverse problem modeling based on PFC simulation exhibits inherent non-uniqueness [12]. Multiple distinct micro-parameter combinations can generate nearly identical macroscopic stress-strain curves. Consequently, a well-fitting inverted parameter set may not be unique. This characteristic is likely a significant factor contributing to the lower \(R^{2}\) observed for \(E_{d}\).