Machine learning-enhanced fully coupled fluid–solid interaction models for proppant dynamics in hydraulic fractures

Statistical analysis of symbolic and CFD-inspired datasets

A detailed statistical analysis was conducted to evaluate and contrast the characteristics of two datasets used to model proppant settling rate (PSR): (i) a symbolic dataset generated using closed-form physical equations, and (ii) a CFD-inspired dataset derived from literature-informed sampling of multiphysics variables. This comparison aims to elucidate the structural and statistical foundations of each dataset, supporting subsequent model development and generalization strategies.

The descriptive metrics summarized in Tables  1, 2 and visualized in the radar plots in Fig. 2 reveal marked differences in the data distributions. The symbolic dataset, constructed under simplified assumptions from Stokes’ law and related formulations, features constrained variability. Notably, gravity and fluid density are constant, reflecting theoretical isolation of variables. In contrast, proppant density, diameter, and viscosity span moderate ranges, aligning with realistic laboratory-scale domains. The corresponding radar plot exhibits sharp peaks for viscosity and particle size, underlining their influence in classical PSR equations.

On the other hand, the CFD-inspired dataset displays a broader and more balanced statistical profile. With variables such as inlet slurry velocity, outlet pressure, elastic modulus, and fracture geometry included, the dataset captures complex subsurface dynamics. Higher standard deviations across these features point to natural operational heterogeneity and multiparametric coupling common in field-scale simulations. The radar plot for this dataset presents a flatter and more distributed pattern, indicating more equitable feature representation. This contrast highlights a critical tradeoff: the symbolic dataset provides clarity and interpretability suited for theoretical validation and symbolic regression, while the CFD-based dataset enables learning in high-dimensional, realistic regimes. Together, they form a complementary foundation for hybrid modeling, bridging idealized physics and data-driven generalization across hydraulic fracturing environments.

Symbolic validation via numerical ODE simulation

To validate the symbolic force-balance formulation, we numerically solved the governing ordinary differential equation (ODE) representing the vertical motion of a single proppant particle in a viscous fluid. The model accounts for gravitational, buoyant, viscous drag, and micro-strain-induced elastic forces acting on the particle. The transient velocity profile (Fig. 3) exhibits a rapid acceleration phase followed by asymptotic convergence to a steady-state terminal velocity, characteristic of drag-limited settling dynamics.

The simulation was conducted for a spherical proppant particle of diameter \(d_p = 0.0006\) m and density \(\rho _p = 3000 {\mathrm{kg/m}}^{3}\), suspended in a fluid with \(\rho _f = 1000 {\mathrm{kg/m}}^{3}\) and viscosity \(\mu = 0.001\) Pa\(\cdot\)s. The equation of motion solved is shown as in Equation 5:

where \(C_d\) is the drag coefficient (corrected for moderate Reynolds numbers), \(k_\epsilon\) is a strain stiffness constant, and \(\epsilon\) represents imposed elastic strain. The term \(k_\epsilon \epsilon\) models elastic retardation arising from micro-strain interactions in confined fracture environments. Specifically, \(k_\epsilon\) is a strain stiffness parameter (in N/m) representing the resistance to deformation encountered by particles as they settle in narrow fracture apertures or near-wall zones. The strain \(\epsilon\) (dimensionless) quantifies the relative deformation imposed by surrounding media, such as fracture walls or adjacent particles. Together, this linear elastic force mimics energy dissipation and velocity damping in response to local structural confinement, a phenomenon relevant to realistic proppant transport scenarios. This formulation is consistent with prior semi-empirical studies on confined particle-fluid interactions³⁸. The drag coefficient varies with Reynolds number Re via Eq. 6:

The particle reached terminal velocity \(\sim\)0.39 m/s within 0.5 seconds, closely matching the analytical Stokes prediction, thus validating the symbolic formulation in the low-to-moderate Reynolds regime. A minor elastic restoring force (\(k_\epsilon \epsilon\)) was included to simulate deformation-induced drag modulation observed in confined fracture networks.

The numerically modeled terminal velocity shown in Fig. 3 is marginally higher than the textbook Stokes prediction. This discrepancy arises due to three main factors: (1) inclusion of drag corrections valid beyond creeping flow (Stokes) assumptions; (2) finite initial velocities leading to faster inertial response; and (3) parametric modeling of micro-strain effects that modify the effective settling resistance. These inclusions were intentional, enabling the symbolic model to emulate behavior observed in realistic proppant transport studies under transitional Reynolds conditions. This validated symbolic model therefore offers a physics-informed, computationally efficient alternative to full-scale CFD–DEM simulations. Its time-resolved outputs may also serve as a benchmark dataset for training surrogate models or physics-informed neural networks (PINNs) tailored for subsurface particle transport.

Parametric sensitivity to diameter and strain

To understand the sensitivity of proppant settling to geometric and mechanical parameters, parametric analyses were performed by varying particle diameter and formation strain independently. Larger particle diameters resulted in substantially faster terminal settling rates, reflecting the quadratic dependence predicted by Stokes’ law. Conversely, increasing formation strain decreased settling velocities, as the strain-coupled elastic resistance effectively opposed gravitational settling. The observed trends reinforces the role of particle-fluid interaction physics in controlling settling behavior.

To explore how particle geometry and formation mechanics influence settling dynamics, a comprehensive parametric study was conducted using the same governing ODEs but with systematically varied inputs. Figure 4 encapsulates the results of these simulations, highlighting the response of settling rate to changes in proppant diameter and formation strain. In the left panel, increasing particle diameter from 0.5 to 1.0 mm led to a marked increase in terminal settling velocity. This is consistent with classical Stokes law predictions, where velocity scales quadratically with diameter, reflecting increased gravitational force and reduced relative surface area available for fluid resistance.

On the right, the influence of formation strain is visualized. In this context, we apply quasi-static strain magnitudes in the range \(\epsilon\) = 0.01 to 0.10 (dimensionless), which represent geometric deformation rather than time-dependent strain rates. The introduction of a strain-dependent resistance force effectively dampens the particle’s descent, mimicking geomechanical feedback from a deforming rock matrix. As the applied strain increases, the net force acting on the particle decreases, thereby lowering its terminal velocity. This nonlinear effect is particularly relevant in unconventional reservoirs, where formation dilation or compaction can significantly alter local flow fields. By incorporating strain as a tunable, dimensionless input in the symbolic model, we bridge the gap between static mechanical deformation and particulate transport, enabling more physically grounded simulations.

Importantly, the dual-panel analysis supports two major conclusions. First, geometric scaling remains a dominant factor in settling behavior, and second, geomechanical coupling cannot be ignored in dynamically evolving formations. The integration of both parameters into the symbolic model allows for more sophisticated feature design when transitioning to machine learning frameworks. Furthermore, these parametric simulations produce synthetic yet physically grounded datasets suitable for training predictive models, ensuring that the learned representations inherit the underlying physics of proppant transport, including nonlinear dependencies and strain-field interactions.

Effect of pressure drop on proppant transport

To extend the symbolic model’s applicability, we evaluated the influence of pressure gradients, rather than absolute pressure drops, across the fracture system. Pressure-induced forces are known to affect proppant mobility, especially under high-pressure injection. Simulations were conducted with pressure gradients of 0.1, 0.5, and 1.0 bar/m to quantify their effects on settling dynamics.

The observed behavior arises from the superposition of a pressure-driven flow field opposing gravity, as shown in Fig. 5a . In the symbolic formulation, this is captured using a Bernoulli-inspired term modulated by fracture deformation. At lower gradients (0.1 bar/m), the drag is negligible, and the particle settles efficiently. At 1.0 bar/m, increased flow significantly reduces settling velocity, highlighting the nonlinear interplay between pressure-driven flow and drag resistance.

This phenomenon has direct implications for field operations. While high-pressure injection facilitates fracture propagation, it may impede proppant transport by suspending particles or diverting them laterally. The model thus highlights a trade-off between fracture extension and proppant placement efficiency. Managing pressure profiles is essential for ensuring particles reach their intended depth before lateral dispersion occurs.

Furthermore, the inclusion of pressure effects underscores the symbolic model’s flexibility and relevance to operational conditions. It also enriches the synthetic dataset for downstream machine learning tasks, where pressure-sensitive dynamics are often underrepresented due to subsurface observability constraints. By simulating pressure-dependent behaviors, the framework provides informative priors that enhance learning accuracy in data-driven models.

Density-driven settling characteristics

The influence of proppant density on settling behavior was evaluated to validate theoretical expectations derived from symbolic formulations. Simulations were conducted for particles with densities ranging from 2500 to \(3500{\mathrm{kg/m}}^{3}\), with results shown in Fig. 5b . Denser particles consistently exhibited deeper settling trajectories and faster descent rates, reinforcing the linear relationship between net gravitational force and density contrast.

This behavior is rooted in fluid mechanics: the net downward force is governed by the density difference between the proppant and the fluid. As density increases, gravitational acceleration grows, while the buoyant force remains constant. Consequently, higher-density particles accelerate more quickly and overcome drag sooner, reaching higher terminal velocities and settling more deeply within a fixed time horizon. The simulated curves clearly stratify according to density, with 3500 \({\mathrm{kg/m}}^{3}\) proppants reaching the greatest depth.

Operationally, these results suggest that selecting denser proppants can enhance penetration efficiency in deep fractures or under turbulent flow conditions. However, such particles may introduce challenges related to transport logistics, suspension stability, and impact-induced formation damage.

The symbolic model includes a geomechanical coupling factor that accounts for formation strain. This correction slightly reduces the computed terminal velocity to reflect deformation effects, such as fracture closure or matrix expansion, which can impede proppant descent. This coupling adds physical realism and increases the reliability of symbolic simulation outputs.

Beyond validation, these synthetic results play a crucial role in informing machine learning models. By simulating a wide range of material properties, the framework ensures that training datasets adequately reflect real-world variability. In doing so, it enhances the generalization capabilities of predictive models tasked with forecasting proppant placement across diverse field conditions.

Combined influence of diameter and density on settling velocity

A symbolic surface model was constructed to jointly examine the effects of proppant diameter and density on settling velocity. The generated response surface, shown in Fig. 6, reveals a monotonic increase in velocity with both parameters. This confirms the theoretical scaling laws, quadratic with diameter and linear with density, and validates the multivariate feature combinations used in downstream machine learning models.

In the Stokes flow regime, diameter has a dominant effect due to its squared contribution to the gravitational term. As particle size increases, the gravitational force (scaling with volume) grows faster than the opposing drag force (scaling with surface area), leading to a nonlinear rise in settling velocity. Simultaneously, denser proppants experience greater net downward force, further boosting descent rates. The interaction between these variables produces a curved, upward-sloping surface that highlights their synergistic effect.

This symbolic simulation serves multiple strategic purposes. First, it affirms the fidelity of the analytical framework, yielding a physically smooth response surface that reflects established fluid mechanics. Second, each surface point acts as a high-resolution synthetic data sample, offering interpretable priors for training machine learning models. Third, it enables informed exploration of the design space, allowing engineers to assess trade-offs. For example, while increasing diameter can improve settling efficiency, it may also raise the risk of bridging or require larger pumps. Likewise, denser materials improve placement but may increase operational costs or pose risks of formation damage.

Unlike traditional CFD, which is computationally expensive for parameter sweeps, the symbolic method generates this insight rapidly with analytical traceability. These results not only support model validation but also ensure that learned behaviors in ML frameworks are grounded in governing physics, facilitating robust predictions in complex fracture environments.

Symbolic pressure drop models under geometric variations

Symbolic expressions capturing pressure drop as a function of particle diameter, density, fracture height, and inlet velocity were evaluated through systematic simulations. The pressure response surfaces under varying conditions, demonstrating consistency with hydraulic principles. Specifically, pressure drop increases quadratically with proppant diameter and density, decreases with fracture height, and rises with inlet velocity, validating the predictive structure of the symbolic model.

Figure 7a reveals that larger and denser proppants cause significantly higher pressure losses. This is attributed to reduced flow area and increased viscous resistance, as larger particles block fluid pathways, and denser materials intensify closure stress and drag effects. The combined influence of these two variables creates steep pressure gradients, suggesting operational limits when using heavy or oversized proppants, especially in tight formations where energy efficiency is critical. Conversely, Figure 7b highlights the strong inverse relationship between fracture height and pressure drop. This follows the cubic dependence of flow capacity on aperture size, where wider fractures dramatically reduce hydraulic resistance. However, increasing inlet velocity produces a proportional rise in pressure demand due to inertial contributions, confirming classical flow resistance trends.

Optimizing proppant selection and injection conditions involves trade-offs: larger particles offer improved transport stability but incur higher pressure penalties, while increasing fracture aperture enhances flow capacity but may be limited by geomechanical constraints. The symbolic framework enables rapid exploration of these multidimensional interactions without requiring computationally expensive CFD simulations. Ultimately, this analysis confirms the symbolic model’s physical fidelity and its practical relevance for guiding parameter tuning in field operations. The resulting surfaces also serve as high-quality priors for training ML models, embedding physics-based generalization into data-driven workflows for more robust proppant transport predictions.

Settling rate and drag force analysis

The settling characteristics of proppants in fracturing fluids are driven by the competing influences of gravity, buoyancy, and drag. To assess how these forces vary with particle properties, a symbolic parametric analysis was conducted across a range of densities and diameters. The results, illustrated in Fig. 8, align well with theoretical expectations derived from Stokesian dynamics.

Figure 8a reveals that the settling rate increases with both proppant density and diameter. This is consistent with the underlying physics, where settling velocity scales linearly with density difference and quadratically with diameter. The smooth gradient in the contour map confirms that denser and larger proppants descend more rapidly due to stronger net gravitational forces and reduced viscous resistance. These insights are particularly relevant for rapid fracture treatments, where minimizing settling time enhances proppant placement efficiency. In contrast, Fig. 8b highlights the drag force behavior, showing that it rises steadily with increasing diameter. This trend reflects the proportionality of drag to surface area and the velocity-dependent nature of fluid resistance. While denser particles settle faster, thereby increasing drag indirectly, the contour plot shows that geometry exerts the dominant control, with drag being more sensitive to diameter than to density.

These trends validate the symbolic model’s predictive accuracy and reinforce its operational value. By mapping settling and drag behaviors against physical properties, the framework offers engineers a practical tool for selecting proppant types based on performance goals. This approach reduces reliance on expensive laboratory tests or full-scale CFD simulations, making it highly suitable for early-stage design workflows. Moreover, the symbolic maps serve as interpretable, physics-grounded features for training data-driven models, improving generalization across variable subsurface conditions.

Predictive performance of symbolic-ML and CFD-ML models

The evaluation of the machine learning models built on symbolic physics-derived and CFD-inspired datasets reveals high predictive performance across all configurations. Figure 9 displays a comparative visualization of model accuracy, showing strong agreement between actual and predicted proppant settling rates (PSR).

The symbolic model, based on first-principles physics, achieved an R\(^2\) score of 0.9934 with a low RMSE of 0.0436 and MAE of 0.0252 m/s, indicating strong generalization in structured, predictable flow regimes. In contrast, the CFD-based model demonstrated a slightly higher R\(^2\) of 0.9941 but with a larger RMSE of 0.2033 and MAE of 0.1603 m/s, reflecting its exposure to more complex, potentially noisier input variations.

The combined stacked ensemble model, which integrates features learned from both symbolic and CFD-inspired datasets, achieved the highest performance overall, boasting an R\(^2\) of 0.9970, RMSE of 0.1801, and MAE of 0.1190 m/s. This result validates the hypothesis that hybridizing physics-based and empirical representations enhances both generalization and robustness.

Across all models as presented in Table 6, the closeness of predicted PSR values to the ground truth confirms the successful capture of nonlinear multivariate dependencies critical to hydraulic fracturing. The symbolic model’s narrower prediction spread confirms that physics-informed feature engineering suppresses noise, while the CFD model’s flexibility showcases its ability to adapt to second-order effects. Together, the stacked model offers a powerful, accurate, and interpretable predictive tool for simulating proppant dynamics under diverse fracturing conditions.

Residual distribution and model robustness

Analyzing the residual distributions provides deeper insight into the generalization capability and error characteristics of the machine learning models. Figure 10a shows the distribution of prediction errors for the combined model. The residuals are heavily concentrated around zero, indicating a high predictive accuracy and confirming that the stacked ensemble model is well-calibrated. The slight right-skewness observed, characterized by a long tail towards positive residuals, suggests occasional underpredictions in settling rate but with an overall low magnitude of deviation.

Furthermore, the residual plot against predicted values, shown in Fig. 10b , demonstrates a random scatter around the zero line without any clear patterns. This random distribution indicates that the model does not suffer from systematic biases across the PSR range, affirming both homoscedasticity and model stability. Notably, there is no visible heteroscedasticity, meaning the error variance remains relatively constant across small and large predicted PSR values. This property is particularly important in practical applications, as it implies that model confidence is consistent regardless of operating conditions.

In the context of hydraulic fracturing, such robustness is crucial. Errors concentrated in specific regimes could mislead operational decisions, such as choosing inappropriate proppant sizes or adjusting slurry rates suboptimally. Instead, the results here suggest that the hybrid symbolic-ML model maintains its reliability across a wide spectrum of flow conditions, supporting its deployment for both low-flow and high-flow fracture environments. These findings also highlight the benefit of integrating physics-derived features, which act as strong priors that guide the model’s learning and suppress erratic behavior even in data-driven layers.

Interpretability through feature importance and PDP analysis

Feature importance and partial dependence analysis offer key interpretive insights into the internal decision-making processes of the machine learning models. Figure 11 illustrates the ranking of input variables based on their contribution to predictive accuracy. The inlet slurry velocity emerges as the most dominant feature across both random forest-based feature importance and permutation tests on the ensemble model. This observation aligns with physical intuition: the inlet velocity directly governs the drag forces and flow regime transitions that critically affect proppant settling.

Secondary influences are observed from fluid viscosity and proppant density, albeit with much smaller contributions. Proppant diameter, surprisingly, plays a relatively minor role, suggesting that within the investigated size range, diameter effects are less nonlinear compared to velocity-driven forces. This reinforces the understanding that in field operations, adjusting slurry injection rates may have a more immediate impact on proppant transport efficiency than altering proppant material properties.

Partial dependence plots (PDPs) in Fig. 12 further elucidate the nature of these relationships. A near-linear positive trend exists for both proppant density and diameter, affirming that larger and heavier particles settle faster. Inlet slurry velocity, however, exhibits a strong nonlinear effect: the PSR remains relatively stable up to about 0.35 m/s, beyond which it rises sharply. This nonlinear transition likely captures a fluid dynamic regime change, potentially from laminar-dominated to inertial or transitional flow, consistent with fracture flow physics.

Together, these findings validate the physical realism of the ML models and demonstrate their interpretability. They affirm that physics-informed feature engineering not only improves predictive power but also embeds operationally meaningful sensitivities that practitioners can leverage for more efficient and optimized fracturing treatments.

Uncertainty quantification of model predictions

Uncertainty quantification is essential for assessing the reliability of machine learning predictions, particularly when models are deployed in decision-critical applications like hydraulic fracturing. In this study, we employed ensemble variance methods to estimate prediction uncertainty, providing a measure of model confidence for each proppant settling rate (PSR) prediction. Figure 13 illustrates the uncertainty bounds using \(2\sigma\) standard deviations around the predicted PSR values.

The analysis reveals that uncertainty tends to increase slightly at higher predicted PSR values. This observation suggests that, while the model performs robustly across most of the PSR spectrum, predictions involving very fast-settling proppants or extreme flow conditions carry more epistemic uncertainty. Such behavior is expected because extrapolating beyond the core training distribution, particularly at larger particle diameters or extreme slurry velocities, naturally incurs greater variability.

Importantly, the majority of predictions fall within narrow uncertainty bounds, indicating that the ensemble model is both well-trained and exhibits low predictive variance across most realistic operating regimes. Figure 14a further confirms that higher model uncertainty correlates slightly with increased RMSE, reinforcing the validity of the uncertainty estimates.

From a practical standpoint, these uncertainty estimates are valuable for risk management in proppant placement optimization. Field engineers can use them to assess the confidence levels of predictive outputs and prioritize intervention strategies when operating near uncertain regions. Moreover, integrating uncertainty-aware models enables probabilistic decision-making, allowing for the incorporation of safety margins in hydraulic fracturing designs. In future extensions, combining this uncertainty quantification with active learning could further enhance model refinement, focusing additional CFD or experimental simulations in regions where model uncertainty remains significant. Thus, uncertainty modeling not only enhances model interpretability but also strengthens the overall trustworthiness of machine learning-assisted hydraulic fracturing workflows.

Feature correlation assessment

Understanding the interdependence among input features is crucial for both model development and physical interpretation. Figure 14b presents the Pearson correlation matrix for all input features considered in the symbolic physics-informed dataset. The matrix reveals that most features exhibit low to negligible pairwise correlations, with the majority of correlation coefficients falling below |0.2|. This independence is a desirable property for machine learning, as it minimizes multicollinearity and helps prevent overfitting during model training.

One moderately positive correlation (\(r = 0.44\)) is observed between proppant diameter and fluid viscosity. This relationship arises because, within the synthetic dataset, variations in viscosity were slightly correlated with variations in proppant size to simulate certain operational scenarios encountered in practice. However, this correlation remains modest and does not significantly impact model generalization.

The orthogonality of features implies that the model can independently learn the distinct effects of each physical parameter on the proppant settling behavior without being biased by strong feature entanglement. For instance, it ensures that the influence of inlet velocity on PSR is learned purely through its impact on drag and acceleration dynamics, rather than being artificially amplified by covariance with other features like proppant density.

From a physics standpoint, the low correlations also validate the symbolic data generation approach, which sampled parameter ranges independently to reflect realistic but uncorrelated experimental conditions. This sampling strategy mirrors controlled laboratory or computational fluid dynamics (CFD) simulations, where variables are independently varied to understand their isolated effects.

The correlation matrix confirms that the feature space is well-structured for machine learning purposes and supports the interpretability of feature importance analyses conducted subsequently. Future work could explore non-linear correlations or mutual information measures to capture potential higher-order dependencies missed by Pearson correlation, thereby enriching the model’s understanding of complex multi-physics interactions.

Validation against literature results

To rigorously evaluate the performance of our symbolic–machine learning (ML) framework, we conducted a comparative analysis against state-of-the-art ML-based approaches reported in peer-reviewed literature for proppant settling rate (PSR) prediction. Table 7 summarizes recent studies across neural networks, support vector regression, recurrent models, and ensemble methods. Accuracy metrics include root mean squared error (RMSE) and coefficient of determination (R\(^2\)), with emphasis on data fidelity and physical realism.

Zhu et al.²² employed a feedforward artificial neural network (ANN) trained purely on geometrical and fluid parameters, achieving an R\(^2\) of 0.972 and RMSE of 0.103. While their model demonstrated reasonable accuracy, it lacked incorporation of physical constraints, thereby reducing generalizability. Mirvakili et al.⁴⁹ extended this approach to non-spherical particles in Newtonian and power-law fluids, reporting RMSE = 0.048 using 361 samples.

More recent work by Lv et al.⁵⁰ introduced shape and orientation features into ANN modeling, obtaining RMSE = 0.026 and 8.6% average relative error, albeit still relying heavily on dense training data. Meanwhile, Tang et al.⁵¹ used a multifidelity neural network (MFNN) to fuse low- and high-fidelity datasets, reducing RMSE from 0.30 to 0.2, a valuable approach under limited experimental availability.

For time-series field data, Hou et al.⁵² applied gated recurrent units (GRU) to predict dynamic proppant pressure (DPP), achieving an RMSE of 7.66 after physics-aware correction. Sharma et al.⁵³ investigated pressure drop prediction in multiphase slurry transport using Random Forests and Neural Networks, reporting R\(^2\) = 0.964, and highlighting the importance of flow rate and mixture density via sensitivity analysis.

Compared to these methods, our symbolic–ensemble approach in Table 8, achieves RMSE = 0.046 and R\(^2\) = 0.997, representing state-of-the-art performance. The integration of physically interpretable symbolic features (from drag, buoyancy, and strain) with ensemble generalization contributes to accuracy, robustness, and transferability across unseen parametric regimes. This affirms the framework’s utility as a scalable, physics-grounded alternative to costly CFD–DEM simulations or black-box ML models.

Study’s practical implication & limitations

The hybrid symbolic–ML framework developed in this study provides a computationally efficient alternative to full-scale CFD–DEM simulations for predicting proppant settling dynamics in hydraulic fractures. By integrating first-principles symbolic models with ensemble learning, the approach retains physical interpretability while achieving high predictive accuracy (R\(^2\) = 0.997 for the stacked model). This makes the framework particularly useful for rapid screening of proppant types, fluid viscosities, and fracture geometries in field-scale design scenarios. Additionally, the symbolic derivation of settling and drag laws enables quick parametric exploration without requiring extensive meshing or fluid–solid interaction solvers.

Despite its strengths, the approach presents several limitations. First, the CFD-informed dataset was generated synthetically and benchmarked against known physical trends but not validated against real experimental data. Future work should include laboratory- or pilot-scale flow loop validation to verify the generalizability of the model to complex field conditions. Second, the symbolic model assumes Newtonian flow, spherical particles, and uniform pressure gradients, which may not fully capture the subtleties of shear-thinning fluids, particle agglomeration, or stratified fracture zones.

Moreover, although ensemble ML models offer robustness, they inherently lack temporal resolution and mechanistic coupling unless explicitly modeled. The current framework is predictive in nature but not prescriptive; it forecasts settling rates based on given inputs but does not simulate transient dynamics or back-calculate optimal input conditions. Finally, while the symbolic model supports interpretability, the stacked ensemble can become opaque without further use of explainable AI tools such as SHapley Additive exPlanations (SHAP) or Local Interpretable Model-agnostic Explanations (LIME).

In all, this study presents a rudimentary step toward scaling physics-informed proppant transport modeling while highlighting the importance of integrating empirical validation and multi-physics generalization in future extensions.