Advanced Mathematical Modelling of Leaching Processes in Porous Media: An Averaging Approach

Leaching, a method characterized by the selective dissolution of minerals from rock through chemical reactions between a leach solution and ore, has gained prominence in the past quarter-century, noted for its cost-effectiveness and minimal environmental footprint. This technique, predominantly adopted in Eastern Europe, Central Asia, and the USA, plays a pivotal role in the extraction of valuable substances. The relevance of mathematical modeling in the leaching process lies in its capacity to precisely predict process outcomes and optimize parameters, thereby enhancing efficiency and cost-effectiveness in substance extraction. These models provide insightful comprehension of the underlying chemical and physical phenomena, which is indispensable for the evolution of current technologies and the development of new ones.

Averaged mathematical models, particularly in the context of rare metal leaching, offer a pragmatic approach for process optimization and control. By focusing on average parameters, these models facilitate predictions of general trends and outcomes, a crucial aspect for decision-making in production settings. In scenarios involving complex chemical and physical processes, such models reduce computational demands and simplify analyses, without significantly compromising the integrity of the essential characteristics of the process. Despite their utility, it is crucial to acknowledge the limitations of averaged models in capturing the intricacies and heterogeneities of the leaching process, which may affect the precision of predictions in certain contexts.

In the sphere of deformable porous media with double porosity, the application of averaged models has been instrumental [1]. These models adeptly segregate various mechanical-hydrodynamic effects, including peristaltic and shear deformation-induced flow. The introduction of new stresses in the deformation equation, attributable to cross effects and matrix relaxation, exemplifies the nuanced understanding afforded by such models. Their primary merit lies in the detailed portrayal of interactions within highly heterogeneous porous media, encompassing both mechanical and hydrodynamic factors. Consequently, averaged models emerge as formidable tools for analyzing and simulating complex processes in deformable porous media, balancing accuracy and practical applicability. Continuing the exploration of averaged mathematical models in the context of leaching processes, Panfilov et al. [2] focused on numerical simulation within a multiscale fractured-porous medium. The essence of this model lies in the averaging of complex media parameters, such as porosity and hydraulic conductivity, thereby facilitating a more streamlined and efficient mathematical representation. This model is distinctive for its portrayal of single-phase flow in a weakly compressible fluid, employing integro-differential equations to compute the average pressure across various sections of the medium.

An averaged model gains significance in the realm of highly dynamic processes. This model finds its application in scenarios like underground gas storage facilities, where alternating phases of gas injection and withdrawal occur. Such environments demand a model capable of adapting to rapid changes, making this averaged approach particularly relevant [3]. Bués and Panfilov [4] adopted a mathematical method to address the leaching of minerals. A differential equation with a delay element encapsulates the transportation of a dissolved substance through a porous medium. This approach is noted for its capacity to provide analytically robust solutions, effectively capturing the systems where current states are influenced by preceding conditions. In scenarios involving delayed decision impact, such as in leaching processes, this method offers a nuanced understanding of temporal dependencies. In-situ leaching (ISL) demonstrates clear advantages over traditional mining methods for non-ferrous metals. These advantages include mitigated changes in underground areas, reduced radiological contamination risk, optimized production delay time, and diminished energy consumption. This method underscores the progressive shift towards more sustainable and efficient mining practices [5].

The numerical analysis presented by Mukhanov et al. [6] delves into the leaching of non-ferrous metals, concluding that the overall resistance of the process is a combination of diffusion and chemical action resistances. It is highlighted that when one stage exhibits significantly higher resistance than others, it becomes the limiting factor for the process's progression. For instance, in cases where a solid product forms a dense layer on the surface of a dissolved mineral, such as uranium, diffusion through this layer becomes the rate-limiting stage. Under these circumstances, the influence of leaching conditions on the rate is governed by the principles of internal diffusion. In the realm of ISL for uranium mining, mathematical modeling serves as an essential tool for accurately delineating the intricate physicochemical processes, including reaction dynamics and hydrodynamic transport [7]. Despite its efficacy in replicating the nuances of leaching, this method necessitates simplifications and case-specific settings, rendering its practical application challenging in the absence of supplementary experimental data.

The research conducted by Simon et al. [8] underscores the need to consider various phenomena in numerical modeling, such as geochemical reactions, reaction kinetics, and the rate of hydrodynamic transfer relative to reaction kinetics. In this study, experiments were conducted to simulate leaching solutions and calibrate the pathways of geochemical reactions, along with their kinetic laws. This research provides insights into processes like colmatage (clogging due to sedimentation) and suffusion (erosion of fine material). Furthermore, a model was developed by Mints [9] to ascertain the concentration of solid suspensions in a fluid and the saturation density of a porous medium. This model facilitates a deeper understanding of the dynamics within the medium. The numerical modeling outcomes of several experimental studies provide a comprehensive view of the subject [10, 11].

Additionally, Khuzhaerov [12] addressed the scenarios where pores were not completely obstructed by the sedimentation of suspension particles. In the context of this model, the filtration of a multicomponent suspension is examined, focusing on the effects of colmatation and suffusion on each component, in relation to the concentration of other components. This analysis yields valuable insights into the interplay of various factors within the leaching process. The innovative approach proposed by De Silva and Ranjith [13] explores new methods of fracture stimulation. By injecting a modified silent crack mitigation agent (SCDA) into the target rock through an injection well, this study demonstrates an effective optimization technique. The simulation results reveal that multi-stage charging amplifies the fracture density of the target rock, thereby enhancing its permeability for fluid flow without necessitating additional injection wells. The research conducted by Van der Lee [14] highlights a crucial prerequisite for the application of mathematical modeling in the leaching process of non-ferrous metals: the permeability of these metals to the leaching fluid. This aspect is fundamental to the effectiveness of any mathematical model in simulating leaching processes. A numerical implementation of algorithms was presented by Zhumali and Reshetova [15], focusing on underground leaching in a two-dimensional framework. This study explores the positioning of soil leaching boundaries under varying parameter values within the differential equation system, alongside examining acid concentration at different temporal points. Such an analysis is vital for understanding the spatial and temporal dynamics of the leaching process.

Mukhambetzhanov et al. [16] introduced a mathematical model that integrates the effects of capillary pressure and gravitational forces, constructed on a balanced monotonic finite-difference scheme. This approach is instrumental in capturing the nuanced interactions within the leaching process. Similarly, the effective numerical modeling proposed by Escandon-Panchana et al. [17] facilitates the accurate prediction of oil production productivity, with potential applicability to mineral extraction in porous media. These models underscore the adaptability and precision of mathematical modeling in diverse applications. Applied studies, such as the study conducted by Cohen et al. [18], demonstrate the practical utility of mathematical modeling in achieving desired production outcomes. Depending on the specific context, researchers can leverage both averaged and detailed models to yield a more accurate analysis of rare metal leaching processes. This flexibility in model selection allows for tailored approaches to suit various research and industrial requirements. Moreover, El Haroui et al. [19] delved into the thermal solution mixed convective flow around a vertical wall in a porous medium, elucidating the impact of free convection, medium permeability, and thermal conditions on heat and mass transfer. Such research underscores the significance of modeling diverse physical processes in porous media, contributing to a comprehensive understanding of these complex systems.

The overarching aim of this study is to develop an improved mathematical model to optimize the leaching process of traditional metals, thereby enhancing efficiency and accuracy in both industrial production and academic research. The application of this model, grounded in averaging methods, facilitates the analysis and simplification of complex chemical and physical processes in porous media. By preserving the main characteristics of these processes, the model ensures sufficient accuracy for critical production decision-making.

The formulated assumptions, coupled with a straightforward mathematical fact proposed by Zhikov et al. [20], underpin the derivation of the averaged dynamic equation of motion for pressure $p(x, t)$. Utilizing the representations (7) and (3) for solutions $p^{\varepsilon}(x, t)$ and $\chi^{\varepsilon}(x, t)$, the equation is written in the form of Eq. (7).

$\int_{\Omega_T}\left(-\alpha_T \mathrm{~V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \cdot \frac{\partial \psi}{\partial t}+\chi\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)(\nabla p(\mathrm{x}, t)+\right.\left.\nabla_y P\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \cdot \psi\right) d x d t=\mathrm{o}(\varepsilon),$                (12)

This dynamic equation, integral in mathematical modeling, delineates changes in parameters such as pressure, temperature, substance concentration, or fluid flow velocity within the pores and cracks of porous materials over time. These equations are pivotal in describing processes like diffusion filtration or chemical reactions in porous media, encompassing soils and geological formations.

For the purpose of analysis, functions of a specific form $\psi=\psi(x, t)$ are chosen as test functions. When approaching the limit for $\varepsilon \rightarrow 0$, the identity expressed in Eq. (13) is derived.

$\int_{\Omega_T}\left(-\alpha_T \mathrm{~V}(\mathrm{x}, t) \cdot \frac{\partial \boldsymbol{\psi}}{\partial t}+(m(\mathrm{x}, t) \nabla p(\mathrm{x}, t)+\right.\left.\left.\int_Y \nabla_y P(\mathrm{x}, \mathrm{y}, t) d y\right) \cdot \psi\right) d x d t=0,$                (13)

where, $v(x, t)=\int_Y V(x, y, t) d y, m(x, t) \int_Y \chi(x, y, t) d y$. The quantity $v(x, t)$ is appropriately termed as the velocity of the continuous medium, and $m(x, t)$ as the porosity of the continuous medium.

The method involves discarding derivatives of the test function through integration by parts. Additionally, leveraging the fact that the integral identity $\int_{\Omega_T} \varphi\left(\frac{\partial a}{\partial t}+\nabla \cdot u\right) d x d t=0$ for arbitrary functions $\varphi$ equates to the equation $\frac{\partial a}{\partial t}+\nabla \cdot u=$ 0 , the study progresses to the differential equation shown in Eq. (14).

$\alpha_T \frac{\partial \mathrm{v}}{\partial t}=-\mathbb{A} \cdot \nabla p$,                (14)

where, $\mathbb{A}(x, t)=m(x, t) \mathbb{I}+\int_Y \nabla_y P(x, y, t) d y$.

This phase of the analysis involves the reintegration of the integral identity. For this purpose, functions of the form $\psi=$ $\varepsilon \nabla\left(\psi_0(x, t) \psi_1\left(\frac{x}{\varepsilon}\right)\right)$, where $\psi_1(y)$ is a smooth 1-periodic function in variable $y$. The progression to the limit for $\varepsilon \rightarrow 0$ yields the identity expressed in Eq. (15).

$\int_{\Omega_T} \psi_0(\mathrm{x}, t)\left(\int_Y \chi(\mathrm{x}, \mathrm{y}, t)\left(\nabla p+\nabla_y P\right)\right..\left.\nabla_y \psi_1(\mathrm{y}) d y\right) d x d t=0$                (15)

The initial stage of reintegration, conducted over variables $(x, t)$, results in the identity $\int_Y \chi(x, y, t)\left(\nabla p+\nabla_y P\right)$. $\nabla_y \psi_1 d y=0$. Subsequently, a second reintegration over variable $y$ leads to the derivation of the differential equation $\nabla \cdot \chi\left(\nabla p+\nabla_y P\right)=0, y \in Y$, interpreted in the context of distribution theory. This equation is conceptualized as a periodic boundary value problem, featuring periodicity conditions for both the solution and its derivatives at the boundaries of the region, which coincide with the boundary $Y_f$ of the unit cube $Y$.

$\Delta_y P=0, \mathrm{y} \in Y_f(\mathrm{x}, t),$                (16)

$\left(\nabla p+\nabla_y P\right) \cdot \mathrm{n}=0, \mathrm{y} \in \gamma(\mathrm{x}, t),$                (17)

where, $Y_f(x, t)=\{y \in Y: \chi(x, y, t)=1\}$, and $n$ is the unit normal vector to the unknown boundary $\gamma(x, t)$ in domain $Y_f(x, t)$.

However, the boundary value problem defined by Eqs. (16) and (17) is not a closed system, as the boundary $\gamma(x, t)$ is subject to temporal and spatial variations and requires determination. Assuming temporarily that the area $Y_f(x, t)$ is known for all $(x, t) \in \Omega_T$, the solution $P(x, y, t)$ to this problem is represented by Eq. (18).

$P(\mathrm{x}, \mathrm{y}, t)=\sum_{i=1}^3 P^{(i)}(\mathrm{x}, \mathrm{y}, t) \frac{\partial p}{\partial x_i}==\left(\sum_{i=1}^3 P^{(i)}(\mathrm{x}, \mathrm{y}, t) \mathrm{e}_i\right) \cdot \nabla p,$          (18)

This formula integrates the standard basis $\left(e_1, e_2, e_3\right)$ of a rectangular coordinate system with solutions $P^{(i)}, i=1,2,3 \ldots$ to the periodic boundary value problem detailed in Eqs. (19) and (20).

$\Delta_y P^{(i)}=0, \mathrm{y} \in Y_f(\mathrm{x}, t),$          (19)

$\left(\mathrm{e}_i+\nabla_y P^{(i)}\right) \cdot \mathrm{n}=0, \mathrm{y} \in \gamma(\mathrm{x}, t)$.          (20)

Considering a unit tensor, denoted as $\mathbb{I}$, and a second-rank tensor $a \otimes b$ operating according to the rule $(a \otimes b) \cdot c=$ $a(b \cdot c)$, the following relation is established as shown in Eq. (21).

$\mathbb{A}(\mathrm{x}, t)=\left(m \mathbb{I}+\int_Y \chi \sum_{i=1}^3 \nabla_y P^{(i)} \otimes \mathrm{e}_i d y\right)$          (21)

The derivation of the remaining averaged equations follows a similar methodology. Incorporating representations Eqs. (7)-(10), the form changes as indicated in Eq. (22).

$\int_{\Omega_T}\left(\chi\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)\left(\rho_f-\rho_s\right) \frac{\partial \varphi}{\partial t}\right.\left.+\rho_f \mathrm{~V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \cdot \nabla \varphi\right) d x d t=\mathrm{o}(\varepsilon).$          (22)

Functions of a specific form $\varphi=\varphi(x, t)$ are chosen as test functions, and upon progressing to the limit for $\varepsilon \rightarrow 0$ and subsequent reintegration, the result is captured in Eq. (23).

$\nabla \cdot \mathrm{v}=\frac{\left(\rho_s-\rho_f\right)}{\rho_f} \frac{\partial m}{\partial t}\cdot$          (23)

For $\varphi=\varepsilon \varphi_0(x, t) \varphi_1\left(\frac{x}{\varepsilon}\right)$ from Eq. (22), the microscopic continuity equation, as represented in Eq. (24), is deduced.

$\nabla_y \cdot \mathrm{V}=0, \mathrm{y} \in Y$.          (24)

This equation is of paramount importance as it delineates the conservation of mass or energy at the microscopic level, accounting for the behavior of individual particles or elements within small volumes.

To elucidate the macroscopic behavior of the reagent concentration, representations Eqs. (7)-(9) are employed, leading to the formation of the integral identity as depicted in Eq. (25). The concentration of the reactant is a crucial parameter, governing the distribution and quantity of chemicals in the pores. This aspect is instrumental in understanding and controlling various chemical processes like filtration, adsorption, and catalytic reactions within porous media.

$\begin{gathered}\int_{\Omega_T} \chi\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)\left((c(\mathrm{x}, t)+\delta) \frac{\partial \xi}{\partial t}-\left(-c(\mathrm{x}, t) \mathrm{V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)+\right.\right. \\ \left.\left.\alpha_c\left(\nabla c(\mathrm{x}, t)+\nabla_y C\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)\right) \cdot \nabla \xi\right)\right) d x d t+ \\ \int_{\Omega} \chi_0(\mathrm{x})\left(c_0(\mathrm{x})+\delta\right) \xi(\mathrm{x}, 0) d x=\mathrm{o}(\varepsilon) .\end{gathered}$         (25)

Similar to the approach taken with the dynamic equation, the diffusion-convection equation and its initial condition are articulated in Eqs. (26) and (27).

$m \frac{\partial c}{\partial t}+\mathrm{v} \cdot \nabla c=\alpha_c \nabla \cdot\left(m \nabla c+\int_Y \chi \nabla_y C d y\right)-\left(\frac{\varrho_s}{\varrho_f} c+\delta\right) \frac{\partial m}{\partial t}$,                 (26)

$c(\mathrm{x}, 0)=c_0(\mathrm{x})$,                 (27)

The diffusion-convection equation encapsulates the process of matter transfer, integrating the mechanisms of diffusion (the random distribution of particles) with convection (their movement under the influence of a flow).

The progression to the limit in the same identity, utilizing test functions of a specified form $\xi=\varepsilon \xi_0(x, t) \xi_1\left(\frac{x}{\varepsilon}\right)$, culminates in the microscopic equation of diffusionconvection $\nabla \chi \cdot\left(\nabla c+\nabla_y C\right)=0, y \in Y$.

When examining the boundary value problem defined by Eqs. (16) and (17), it becomes evident, as denoted in Eq. (28), that certain conditions apply.

$C(\mathrm{x}, \mathrm{y}, t)=\sum_{i=1}^3 P^{(i)}(\mathrm{y}) \frac{\partial c}{\partial x_i}=\left(\sum_{i=1}^3 P^{(i)}(\mathrm{y}) \otimes \mathrm{e}_i\right) \cdot \nabla c$,                (28)

Consequently, Eq. (26) evolves into the finalized form of the diffusion-convection equation for the concentration of the reagent, as expressed in Eq. (29).

$m \frac{\partial c}{\partial t}+\mathrm{v} \cdot \nabla c-\nabla \cdot\left(\alpha_c \mathbb{A} \cdot \nabla c\right)=-\left(\frac{\varrho_s}{\varrho_f} c+\delta\right) \frac{\partial m}{\partial t}\cdot$                (29)

In a parallel vein, the equations of diffusion-convection involving solid skeleton matter are deduced, following a similar methodology, and are depicted in Eq. (30).

$m \frac{\partial c_i}{\partial t}+\mathrm{v} \cdot \nabla c_i=\frac{\rho_s}{\rho_f}\left(\gamma_i c^{v_i}-c_i\right) \frac{\partial m}{\partial t}\cdot$                (30)

The integral identity corresponding to representations Eqs. (7)-(9), as indicated in Eq. (31), reveals the most significant challenge in this context is the transition to the limit for $\varepsilon \rightarrow 0$, particularly concerning the term

$\begin{gathered}I_{\varepsilon}=\int_{\Omega_T} \psi c^{v_i-1} \mathrm{~V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \cdot \nabla_y C\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) d x d t \cdot \int_{\Omega_T} \chi\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)\left(\left(c_i \frac{\partial \psi}{\partial t}\right.\right. \\ \left.+\left(c_i-\frac{\rho_s \gamma_i}{\left(\rho_s-\rho_f\right)} c^{v_i}\right) \mathrm{V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \cdot \nabla \psi\right)-\frac{\rho_s \gamma_i v^i}{\left(\rho_s-\rho_f\right)} \psi c^{v_i-1} \mathrm{~V}\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right) \\ \left.\cdot\left(\nabla c+\nabla_y C\left(\mathrm{x}, \frac{\mathrm{x}}{\varepsilon}, t\right)\right)\right) d x d t=\mathrm{o}(\varepsilon), \\ i=1, \ldots, n,\end{gathered}$                (31)

Given that for $\psi=\psi(x, t)$, by virtue of Eq. (24), it is established that $\lim _{\varepsilon \rightarrow 0} I_{\varepsilon}=I=\int \Omega_T \psi c^{v_i-1}\left(\int_Y V \nabla_y C d y\right) d x d t$ and, $\int_Y V \nabla_y C d y=-\int_Y C \nabla_y V d y=0 I=0$. This finding is instrumental in progressing towards the derivation of the final boundary condition on the desired boundary $\gamma(x, t)$ of the domain $Y_f(x, t)$.

The third fundamental assumption, as articulated in Eq. (32), is recalled.

$\frac{\beta}{\delta}=\varepsilon \lambda$.                (32)

Utilizing representation Eq. (2), the function $a$ is chosen to be as defined in Eq. (33), and the specifications are outlined in Eq. (34) and $n$ represents the outer unit normal vector to the boundary $\gamma(x, t)$.

$a^{\varepsilon}=a\left(x, \frac{x}{\varepsilon}, t\right)$,                (33)

$\mathrm{a}(\mathrm{x}, \mathrm{y}, t)=\mathrm{n}(\mathrm{x}, \mathrm{y}, t), \mathrm{y} \in \gamma(\mathrm{x}, t)$               (34)

As $\varepsilon$ approaches zero, this leads to the integral identity shown in Eq. (35).

$\int_0^T \int_{Y_f(\mathrm{x}, t)}\left(\frac{\partial}{\partial t}\left(\left(\zeta_1(\mathrm{y}, t) \cdot \mathrm{a}(\mathrm{x}, \mathrm{y}, t)\right)\right.\right.\left.\left.+\nabla_y \cdot\left(\lambda c^v(\mathrm{x}, t) \zeta_1(\mathrm{y}, t)\right)\right)\right) d y d t=-\int_{Y_f(\mathrm{x}, 0)} \zeta_1(\mathrm{y}, 0) \cdot \mathrm{a}(\mathrm{x}, \mathrm{y}, 0) d y$.               (35)

The subsequent conversion of the volume integral into a surface integral, followed by reintegration, results in the final boundary condition expressed in Eq. (36). This condition pertains to the value of the boundary $(\gamma)$ velocity $V_v$ in the direction of the normal $v$ to this boundary and the initial condition detailed in Eq. (37). These conditions are pivotal for determining the position of the boundary $\gamma$ for $t>0$.

$V_v(\mathrm{x}, \mathrm{y}, t)=\lambda c^v(\mathrm{x}, t), \mathrm{y} \in \gamma(\mathrm{x}, t)$,               (36)

$\chi(\mathrm{x}, \mathrm{y}, 0)=\chi_0(\mathrm{x}, \mathrm{y})\left(\gamma(\mathrm{x}, 0)=\gamma_0(\mathrm{x})\right)$.               (37)

In the presented model, a small parameter, closely associated with the average dimensionless pore diameter in porous media, is identified as the cornerstone of the methodology. This parameter is instrumental in the process of averaging the system of equations, a critical component of the mathematical framework. Such averaging is pivotal for the generalization and simplification of the original mathematical model. It renders the model more practical for analytical and computational applications while preserving the essential properties and characteristics of the original process. This approach markedly reduces computational complexity, an aspect of paramount importance in dealing with intricate chemical-physical processes like leaching. A reduction in computational demand not only streamlines the calculation process but also enhances the model's applicability in diverse scenarios, particularly those with constrained computational resources.

Furthermore, the model's simplified nature aids in elucidating the core mechanisms of the leaching process. This insight is invaluable for researchers and engineers specializing in hydrodynamic process modeling within porous media. By implementing the averaging method, the model successfully maintains overall accuracy in system behavior predictions, considering the average parameter values. Consequently, the model emerges as a potent tool for forecasting general trends and outcomes in leaching processes, facilitating more efficient planning and management. It also contributes significantly to the development and optimization of new methodologies within this domain.

The results derived from this mathematical modeling underscore the efficacy of utilizing averaged models for analyzing leaching processes. These models adeptly simplify complex chemical-physical processes while retaining their fundamental characteristics, thereby playing a crucial role in accurately predicting leaching processes.