Abstract

Anisotropic and heterogeneous solids, comprising a minimum of two or more elements with different properties, appear pervasively in rock materials including the pore structure and mineral composition of quartz and granite and usually have extensive applications in construction aggregates, dimension stone, geotechnical engineering, and petrology/geology. The elastic stress fields of inhomogeneous materials or composites inevitably change due to the presence of heterogeneous microstructure under applied external conditions, and hence the total mechanical and physical properties are affected by the temperature and pore pressure fluctuations beneath the surface. In this paper, the thermoporoelastic (TPE) approach on the basis of eigenstrain concept is introduced to predict the stress fields in fluid-saturated porous geologic materials like hydrocarbon reservoirs or aquifers by accounting for the coupling between thermal, poroelastic, and mechanical effects. It is an extension of the micromechanical theory that also incorporates thermal effects like pore pressure changes and temperature alternations. In addition, the TPE approach provides an important multiphysical modeling tool for understanding subsurface fluid-rock interactions and stress states in applications like unconventional oil/gas and geothermal energy.

1. Introduction

Eshelby was a pioneer in micromechanics who insightfully introduced the concept of transformation strains to describe the effect of point defects and dislocations on the elastic field. The Eshelby inclusion problem was addressed to solve the elastic fields with respect to the stresses and strains around an ellipsoidally shaped region that undergoes a spontaneous change of form inside a full media. He presented the thought experiment involving imaginary cutting, straining, and welding operations [1] and developed his theory to explore the problem of matrix and inclusion with different elastic constants via the harmonic potential functions [2]. His highly influential works have applications in a spectrum of engineering fields which encompass materials science, continuum mechanics, and geological engineering. His model can be utilized to investigate the deformation of rock systems and minerals due to inherent cracks, pores, and grains. Later, Mura extended Eshelby’s theory on the concept of eigenstrains to study the mechanical behaviors of materials at the microscopic level and handled various problems relating to inclusions, dislocations, cracks, composites, and polycrystals [3]. Biot, who is the founder of the theory of poroelasticity, and Willis discussed how to choose a suitable combination of measurements including shear modulus, jacketed and unjacketed compressibility, coefficient of fluid content, and porosity to determine the coefficients for an isotropic system and established the measurement and interpretation of the elastic coefficients of Biot’s theory of deformation in a porous elastic solid containing a compressible fluid [4]. This paper also contains extension methods to deal with anisotropic systems, linear systems, and nonlinear systems with proper stress definitions. The analytical and numerical solutions are developed to deal with the issues of gas-water production, geothermal system, and low-permeability reservoir in the petroleum industries and related engineering fields [57].

Oil, gas, and geothermal production can induce earthquakes by altering pore pressures, temperatures, and stresses within and around reservoirs. Near reservoir margins or in high-pore pressure gradient areas, dilatant fracturing and normal faulting are always promoted in extensional settings, which can enhance fracture permeability in adjacent tight rocks [8]. Rudnicki generalized the Walsh method of the pore pressure alternation expressed in terms of the transformation strains [9]. The results can be widely applied in computing the stress, strain, and surface displacement fields produced by fluid mass or extraction and pore pressure alternation within a resource beneath the ground for the purposes of hydrocarbon production and underground storage, aquifer management, and carbon sequestration. The interaction mechanism between two ellipsoidal poroelastic inhomogeneities, whose methodology is devised by Moschovidis and Mura [10], was provided to elucidate the stress perturbations around the inhomogeneities with the assistance of the equivalent inclusion method and the higher rank of Eshelby’s tensors [11]. The proposed solution is applied in calculating the stress and strain fields around the arbitrarily orientated hydrocarbon resources caused by fluid injection or withdrawal [12]. An extension problem related to the poroelastic damage to the brittle rock failure was addressed under the framework of microstructural and hydromechanical methods [13]. The randomly distributed spheroidal microcracks with oblate shapes and nanopores were considered within the computational model. This provides a systematic approach to incorporating microstructural mechanisms in brittle fractures. Furthermore, the nonlinear and heterogeneous problem of layered cylindrical inhomogeneity in the field of depleting resource system with concerning the caprock formation was solved through Kelvin’s solution and the discrete collocation fixed point iteration approach [14]. The stress evolution with considering 4D analysis of strain-dependent elasticity and application in subsalt carbonates was investigated in comparison with FEM results. The special geometry of a thin disk-shaped resource subject to depressurization was explored via the thermoporoelastic inclusion model [15]. In addition, the related topics on the stress response to pore pressure fluctuations around resources [16, 17] and the nonprobabilistic model to solve the effect of crack size uncertainty [18] were analytically evaluated. In our previous study [19], the induced stress redistributions of the penny-shaped reservoir were investigated via the classical Eshelby inclusion approach. The solved full-space problem on the specific shape of the resource can be extended to treat the associated half-space geostructural problem.

A continuum damage mechanics model was developed to describe the elastic, plastic, and damage behavior of porous rocks [20]. Microcrack and microvoid nucleation and coalescence were introduced within the fracture mechanics framework, and the developed model can be used for hydraulic fracturing simulations in reservoir rocks by eliminating stress singularities and simulating progressive failure without remeshing. The transitions of failure mode on the rock cutting with considering the rock heterogeneities with respect to the microcracks, intergrain cracks, temperature variations, and confining pressure were investigated via numerical methods [2123]. Furthermore, the well integrity of inclusion-matrix systems with concerning the effects of pressurized cracks which are randomly distributed around a rigid inclusion or inclusion-matrix interface was theoretically solved via the theory of eigenstrains [24]. The suggested solution can be applied to a hollow cylindrical casing-cement sheath-formation rock system for corrosion-related problems. Moreover, the ultrasonic application on the basis of the Eshelby-Cheng effective medium theory for porous vertically transversely isotropic media [25], the geomechanical effect of low-temperature CO2 injection via the coupled thermoporoelastic model [26], and the stress analysis and pore pressure variations of hydrocarbon-bearing formation by the reservoir geomechanical model [27] were theoretically investigated.

The fundamental solution for a continuous line source injecting into a poroelastic reservoir bounded by impermeable elastic layers was derived, which can model subsurface fluid injection or production through a vertical well [28]. The pore pressure field is decoupled and governed by the classical diffusion equation solution for an infinite line source, while the mechanical fields are resolved by an elasticity problem with a body force dependent on the time-varying pore pressure gradient. The Eshelby inclusion theory can be applied to simulate fluid extraction and injection-induced stress perturbations in a porous matrix for half-space-related problems. Employing this theory, approximate analytical solutions for finite-depth resources with rectangular and elliptical geometries under a plane strain condition are derived [29]. The Coulomb failure stress change technique is then utilized to evaluate the fault reactivation potential resulting from the induced stress alternations. Normalized stress change factors and stress variations with different reservoir shapes and depths were determined via the proposed approach. The porothermoelastic problems on the single- and double-inclusions [30] and hydraulic fractures [31] beneath the surface were addressed by using semianalytical and analytical solutions. The fracturing and fault reactivation of the caprock system can also be solved through Eshelby’s theory and the Monte Carlo simulation framework [32]. The closed-form solution for the geometry-simplified reservoir was provided to predict the stress and critical pressure changes due to porosity modulation in conjunction with the Coulomb failure criterion. The related issues of stress concentrations at fault tips and fault length effects of the caprock formation [33] and six strains isolated in the low-permeability layers to examine the influence of the consortium strains [34] were addressed.

2. Fundamental Equations

2.1. Stress State of a Dilatationally Eigenstrained Inclusion beneath the Surface

The materials and methods section should contain sufficient detail so that all procedures can be repeated. It may be divided into headed subsections if several methods are described. A semi-infinitely extended elastic solid, , containing an ellipsoidally shaped inclusion with prescribed porothermal eigenstrain, (simplified as ) is considered. The shear modulus and Poisson’s ratio are denoted as and , and represents Kronecker’s delta. The coefficients of linear thermal expansion and Biot , the corresponding local fluctuations of temperature and pore pressure are involved within the equation. The ellipsoidal domain is expressed in terms of semiaxes , , and and depth location (Figure 1).

The elastic moduli for both matrix and inclusion are denoted as . The stress distributions due to a half-space inclusion can be expressed as follows [35]: where the total strains are The potentials are defined as wherein the term and for the exterior points of subdomain, Ω, and for the interior points when . The corresponding coordinate transformations related to the , , and are detailed in [35]. The terms , , and mean the derivatives with respect to , i.e., ψ/xixj, ϕ/xixj, and .

With the help of the above expressions, the stresses for the interior field are

When the , the stress components for the exterior field may be obtained from Eq. (2).

2.2. Flowchart of Poro-Thermo-Geomechanical Modeling

The presented study outlines a computational approach to model the elastic stress field caused by a porothermo-Eshelby inclusion beneath the surface. The flowchart summarizes the key steps involved in implementing the proposed solution (Figure 2). First, the geometrical parameters defining the shape of the ellipsoidal inclusion are specified. These include the three semiaxes , , and which characterize the size of the ellipsoid. Additionally, the depth location of the inclusion center beneath the surface is input as a model parameter. With the geometry established, the elastic properties of the homogeneous matrix surrounding the inclusion are defined. Young’s modulus, shear modulus, and Poisson’s ratio for the matrix need to be provided. Together, these constants characterize the isotropic elastic behavior of the matrix medium. Next, the porothermal eigenstrains within the inclusion domain account for the changes induced by the temperature variation and pore pressure fluctuation . The linear thermal expansion behavior and the change in fluid content which are related to the interior pore pressure are governed by the thermal expansion coefficient and the Biot coefficient . The porothermal strains are thus obtained by combining the thermoelastic and poroelastic effects.

To find the displacement field produced by the eigenstrained inclusion, a set of elliptic and potential functions must be solved. This is accomplished analytically using explicit equations with proper outward unit normal vectors. The solution provides the three components of the displacement vector at each point in the model domain. From the displacement field, the total strain tensors are calculated at each point by taking the derivatives of the displacement vector. The strain accounts for both the imposed eigenstrains within the inclusion as well as the induced elastic strains in the surrounding matrix. In order to solve the corresponding stress field, the position of the point of interest relative to the inclusion domain must be determined. For interior points, the stress is evaluated directly from Eq. (5) which is derived through the strain and elastic stiffness tensor. For the points outside the inclusion, the stress is obtained from Eq. (6). The stress distribution along any desired line can then be obtained by querying the stress tensor components at a sequence of points along the line. This provides the model output, showing how the porothermal inclusion distorts the stress state along a target line through the domain.

The proposed computational methodology provides an efficient analytical solution to predict the elastic stress field generated by an porothermo-Eshelby inclusion situated within a semi-infinite homogenous, isotropic, linear elastic solid. The technique can be outlined in a nine-step workflow as follows: Step 1 (geometrical definition of inclusion)—the geometry and location of the buried inclusion are defined. The inclusion is modeled as an ellipsoidal domain with an arbitrary shape embedded at a finite depth in the semi-infinite matrix. The major and minor axes of the ellipsoid and its burial depth serve as key input parameters that influence the near-field stresses. Step 2 (property specification of inclusion and matrix)—the mechanical properties of inclusion and matrix material are specified by defining their Poisson’s ratio and Young’s modulus. The homogeneous matrix properties are assumed to be constant. Step 3 (porothermal eigenstrain calculation)—the mismatch strains induced within the inclusion by thermal expansion and pore pressure effects are quantified through an eigenstrain term. This inelastic strain drives the evolution of the elastic stress field. The porothermal eigenstrains are calculated based on inclusion properties. Step 4 (definition of potential functions)—the analytical form of the potential functions is introduced to describe the elastic field induced by the eigenstrain. These potentials satisfy the governing equilibrium equations for the defined problem.

In addition, Step 5 (displacement field solution)—the displacement components throughout the matrix are derived from the potential functions. The displacements quantify the deformations of the inclusion and matrix arising from the buried inclusion. Step 6 (total strain computation)—with the displacement field solution, the total strain tensors are computed by evaluating spatial derivatives. The strain tensors describe the state of localized, multidirectional deformation. Step 7 (interior/exterior point designation)—the spatial location of interest is assessed relative to the domain of the inclusion. Interior points within the ellipsoid and exterior points in the far-field matrix satisfy different equilibrium conditions. Step 8 (stress evaluation at the point of interest)—using the constitutive equations of isotropic linear elasticity, the stress tensors are computed at the designated point based on the interior or exterior field solution. Step 9 (stress distributions output)—by performing steps 7 and 8 at numerous points, the full stress distributions are mapped. The perturbation in the stress field caused by the buried inclusion is visualized over the spatial domain, providing insight into the near-domain thermomechanical effects. In summary, the solution methodology relies on a combination of analytical equations to describe the near-field elastic behavior. By leveraging the computational technique, the stress state induced by a buried thermally and porously expanding inclusion can be modeled efficiently. The step-by-step workflow outlined here provides a roadmap for implementing the proposed approach to gain insight into subsurface thermomechanical processes.

3. The Influential Aspects

The shape of a reservoir, including any alternations made through drilling or production operations, can significantly influence the stress distribution. The location and depth of a reservoir affect the in situ stress state due to regional and local rock loads as well as tectonic forces. Changes to pore pressure from fluid injection/withdrawal will modify the effective stress acting on reservoir rock. Thermal effects from processes like steam injection can also impact stresses. Heating reduces rock strength, causing differential compaction and redistribution of stresses (Figure 3). In summary, geoengineering activities that change reservoir geometry will alter load paths and may cause localized stress concentrations or rotations. Any alternations in these factors may subsequently change stress distributions and orientations in both reservoirs and adjacent formations. The complex interactions between reservoir geometry, location, depth, fluid pressures, and temperatures govern in situ stresses. The thermoporoelastic response of reservoirs and surrounding rocks results in stress changes that must be considered within the geomechanical model. In order to study the stress evolution in engineered reservoirs, the normalized factor is set for the comparative analysis of relative consistency.

3.1. The Geometrical Changes of Reservoir Region

The production process may cause some geometrical changes in the rock and the hydrocarbon region, which may affect the flow behavior and the efficiency of the injection. One of the geometrical changes is the deformation of the rock due to the fluid pressure and stress. The rock may expand, contract, or fracture depending on the properties of the rock and the fluid. The deformation may alter the porosity and permeability of the rock, which are important parameters for fluid flow. The deformation may also change the shape and size of the hydrocarbon region, which may affect the contact area and the displacement efficiency of the injected fluid, as well as the chemical reactivity and stability. Another geometrical change is the phase transition, which may affect the density, viscosity, and compressibility of the hydrocarbons, which are also important parameters for fluid flow. The effect of shape alternation ( varies from 5 to 10) on the stress evolutions of the hydrocarbon region (, ) beneath the ground at the depth of is evaluated (Figures 4 and 5). The coefficients of linear thermal expansion and Biot , the corresponding local fluctuations of temperature and pore pressure are set.

The graphs illustrate three different stress-axis curves, each with varying aspect ratios. The red curve () in Figures 4(a) and 4(b) shows that and display the relatively larger values in the beginning and end than those of the green () and dark yellow curves (). The dark yellow curve has the lowest stress values among the three, starting at a lower point than the red and green curves and increasing more gradually. It is noted that of the dark yellow curve exhibits a relatively low stress in the beginning, gradually decreasing and reaching a bottom at approximately -0.1. The other two curves follow a similar pattern of increasing stress, and their peak is reached at a higher value of around 1.2. As compared to Figures 4(a) and 4(b), normal stress in Figure 4(c) exhibits the same trend at the beginning and end for three kinds of aspect ratios and reaches around 0.1 and a minimum value of -0.9 at the geometrical boundaries for different shapes. Moreover, the shear stress in Figure 4(d) of the red curve starts at a larger point and gradually increases to a peak at approximately 0.025, then decreases to the minimum point at -0.025. The difference in the shapes and slopes of the curves can be attributed to the varying aspect ratios, which affect the distributed stress components along the axis.

In geology, Von Mises stress can be used to analyze the deformation and failure of rocks and soils under the influence of natural or artificial factors such as injection, production, and groundwater flow. By calculating , it can be determined whether the rock and soil mass reach the critical condition of yielding or failure and take corresponding protective measures. The red curve shows a relatively higher stress value of 0.4 at the beginning and then exhibits a more gradual incline as compared to the green and dark yellow curves, which illustrate the stress jumps inside and outside the geometrical boundaries. The dark yellow curve lies between the red and green curves, demonstrating an intermediate value and a moderate incline. The smaller the aspect ratio, the larger the stress jump at the interface between the reservoir and rock formation. The overall trend of the curves shows a progressive increase in stress values as the aspect ratios increase. Understanding these relationships between stresses and aspect ratios can provide valuable insights into the rock system’s behavior during production.

3.2. The Location and Depth Impact

The combination of location and depth determines the initial stress state of the reservoir rock. Local geology can lead to variations where horizontal stresses are dominant even at depth. Accurately characterizing the stress distribution is crucial for managing reservoir productivity and stability during drilling and production operations. The complex interplay between location, depth, and in situ stresses must be understood to effectively engineer reservoir behaviors. In addition, location and depth are two primary controls on the stress profile within a reservoir. Their collective influence arises from regional and local rock loads as well as tectonic forces unique to each area and depth. A careful assessment of these factors provides critical insights into reservoir stress distributions. The effect of various depths (from to ) for the penny-shaped resource (, , and ) on the stress redistributions near the surface is estimated (Figures 6 and 7). The linear thermal expansion and Biot coefficients and fluctuations of local temperature and pore pressure are the same as in Section 3.1 for the purpose of comparations.

The red, green, dark yellow, and blue curves, respectively, represent the stress variations of inclusions located at a depth from to . It can be seen that of the red curve in Figure 6(a), whose inclusion is placed closest to the ground, exhibits a relatively higher value at the beginning and end than those of others, while the stress inside the region shows the greatest difference between the boundary and the center of the inclusion. The same trend is also valid for and distributions of the red curve in Figures 6(b) and 6(c). The blue curve which means the deepest inclusion beneath the surface shows the smoothest transition curve in , , and among the four. The stress components inside the region are not disturbed by the surface effect. As can be seen from Figure 6(d), the shear stress components generally exhibit different paths along with changes in the depth of inclusion. The inclusion located closest to the surface (red curve) experiences a sudden change in shear stress, exhibited as two fluctuations of high and low peak values before and after the boundary. The shear stress component of the deeper inclusion (green curve) illustrates a relatively smoother change as compared to the shallowest one, whereas the dark yellow and blue ones show the same trends with increasing depth.

The graph in Figure 7 displays four different of an inclusion with various depths as compared to Figure 6. The four curves exhibit generally similar trends, where the red curve has slightly higher tensile stress values at the beginning and end compared to the other three curves, while its tensile stress values are lower within the thermoporous region. Notably, the maximum tensile stress values of the four curves appearing at the boundary of the region are approximately the same as 1.9. Within the region, the of the red curve (depth at ) is lowest at the center point of inclusion, while the corresponding values for the other three curves are highest. As the depth increases, the stress values tend to converge, while the inclusion closer to the surface is more heavily influenced by the surface effect. The slopes of the curves differ due to the varying depths of the inclusion, which affect the distribution of around the surrounding rock formation.

3.3. Effect of Thermal Fluctuations

Temperature changes and thermal fluctuations can induce significant stresses in reservoir structures due to thermal expansion and contraction. As the temperature rises, the heated fluid expands, putting outward pressure on the reservoir walls and floors. Rapid cooling of the reservoir has the opposite effect, as contraction of the fluid causes depressurization that changes the stress distribution on structural elements. Repeated cycles of heating and cooling lead to thermal fatigue that can compromise the integrity of the reservoir and surrounding formation over time. Thermally induced stresses also affect reservoirs’ spillways, gates, and outlet controls which must continue to operate reliably despite fluctuating pressures and loads. Therefore, accounting for thermal fluctuations is an important consideration in reservoir management and operations. The response of temperature alternations with varying thermal fluctuations (, -25, 0, 25, and 50) to the stress redistributions is determined in a half-space (Figure 8). The location position of the reservoir is at a depth of , and the other parameters except for the thermal changes are the same as those in Section 3.2 for the comparative analysis.

The graph describes five different stress curves with respect to the red (), green (), dark yellow (), blue (), and purple () of a half-space inclusion produced by thermal fluctuations. In response to the alternations in the interior temperature field, inside and outside the inclusion shows a symmetrical trend of variation, with the maximum values at around 70 occurring at the boundaries. As the temperature changes in the reservoir, the stress jumps at the boundaries are about 30 for each segment of temperature increase in Figure 8(a). Likewise, the curves for exhibit symmetry as well, with the maximum and minimum values differing by approximately 90 in Figure 8(b). The temperature change appears to have a greater impact on the stress redistribution at the reservoir’s boundaries, while inside and outside the reservoir are relatively smaller. Regarding the changes in shear stress in Figure 8(c), as the temperature inside the reservoir continuously increases, its tensile stress value also increases continuously, but near the reservoir boundary, it rapidly decreases to compressive stress. While inside the reservoir, transitions from compressive stress to tensile stress, passing through zero at the center point, and then transitions back from compressive to tensile stress. Similarly, in the case of continuously decreasing reservoir temperature, the changes in are opposite to those of increasing temperature. It is worth noting that the maximum values of occur at the boundary points regardless of the highest temperature increase or decrease. The exhibit the maximum tensile stress of around 95 at and . This graph is important as it helps to understand the effect of thermal fluctuations on the stress levels of a half-space inclusion under different thermal conditions, such as specific heat management properties or analyzing the durability of a reservoir exposed to temperature alternations.

3.4. Pore Pressure-Sensitive Analysis

Before injection, the interior pressure of the reservoir is in its initial state, which depends on the depth and geothermal gradient. During processing, fluid is injected into the reservoir, leading to an increase in pore pressure. This stress release causes the pore throat size to increase and microfractures to open up, enhancing the permeability. In the condition of production, the reservoir pressure is elevated compared to the initial conditions. The increased pore pressure reduces the effective stress acting on the reservoir rock. This stress change causes rock compaction and shifts stress-dependent petrophysical properties like porosity and permeability. The altered pressure regime also impacts the geomechanical stability and changes the propensity for subsidence or induced seismicity. Furthermore, the cyclic changes in pore pressure during fluid flooding lead to reversible alterations of reservoir stress state, pore structure, and petrophysical properties. Proper assessment of these dynamic effects is crucial for optimizing waterflood performance and geomechanical stability. The sensitive analysis of interior pore pressure variations (, -25, 0, 25, and 50) on the stress state of the injected region is determined as compared to the preceding case. Figure 8 depicts a comparative study of the pore pressure effect on the geophysical changes.

The pore pressure alternations varying from to 50 and impacting the exterior and interior stress fields of the resource are explored. As the internal pressure increases, gradually decreases, as shown in Figure 9(a). Due to the influence of the ground, the values at the region center points are relatively smaller compared to those at the region boundaries. However, the trends of are generally the same as the pressure changes in Figure 9(b). The difference is that the jump values at the region boundaries reach the maximum compressive stress when the pressure drops to -50 and reach the minimum compressive stress when the pressure rises to 50. Similarly, the same trend also appears for in Figure 9(c). The difference is that, compared to when the pressure increases, the tensile and compressive states of at the boundaries and slightly outside the domains are slightly higher when the pressure decreases to -50. As can be seen from Figure 9(d), has the maximum values at about 2.1 when the pressure decreases and has the minimum values at about 1.7 when the pressure increases at the interfaces. Due to the difference between the reservoir and rock formation caused by the existence of inherent strains in the inner field and vanishing in the outer field, the maximum also appears at the boundaries.

4. Conclusions

(1)Evaluating the fluid flow, heat transport, interior temperature, and pore pressure fluctuations inside the hydrocarbon energy during production is of great practice. The disturbance effect on the shallow geostructure affecting the geomechanical properties around the rock formation is difficult to solve. It is therefore necessary to adopt the micromechanics method to handle quasistatic and isothermal poroelastic problems(2)Based on Green’s function and coordinate transformations for the semi-infinite isotropic solid, the stress state of a dilatationally eigenstrained inclusion beneath the surface can be derived in terms of potential functions. The approach uses constitutive laws and governing equations that link changes in stress, strain, temperature, and pore fluid pressure of the subsurface reservoir in a porous matrix(3)The thermal and pore pressure fields are coupled through the eigenstrain equation within the domain. The resource located closest to the surface experiences a sudden change in shear stress, which is exhibited as two fluctuations of high and low peak values before and after the boundary. As the depth increases, the stress values tend to converge, while the resource closer to the surface is more heavily influenced by the surface effect(4)TPE models are often complex and computationally intensive, requiring advanced analytical techniques like potential theory for a better understanding of rock failure/fracturing, reservoir compaction, surface subsidence, and induced seismicity that can occur with fluid injection/withdrawal. It is recommended that Eqs. (5) and (6) can be employed to predict the stress changes resulting from operations like hydrocarbon production or injection, geothermal energy extraction, and CO2 sequestration

Data Availability

The processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work is financially supported by the Natural Science Foundation of Chongqing, China (Nos. cstc2021jcyj-msxmX1109 and cstc2020jcyj-msxmX1035) and the Doctoral Program of Chongqing, China (No. CSTB2022BSXM-JCX0167). The authors would also like to acknowledge the support from the Youth Project of Science and Technology Research Program of the Chongqing Education Commission of China (Nos. KJQN202203206 and KJQN202103223).