Fractal Dynamics Modeling of Spatiotemporal Pore Structure Evolution in Tight Reservoirs

1. Introduction

The pore structure of tight reservoirs fundamentally controls hydrocarbon recovery, fluid transport, and subsurface carbon storage, However, its inherent complexity and heterogeneity continue to pose formidable challenges for reservoir characterization [1,2,3]. Shales, coals, and carbonate formations exhibit multi-scale pore networks ranging from nanometer-scale micro-pores to interconnected fractures, with pore heterogeneity directly influencing adsorption capacity, permeability, and reservoir quality. Recent advances in experimental characterization, including low-pressure gas adsorption, nuclear magnetic resonance (NMR), mercury intrusion porosimetry (MIP), and high-resolution imaging, have revealed that pore systems follow fractal scaling laws, offering a quantitative descriptor of their structural irregularities [4,5,6,7,8,9].

Fractal geometry has thus become a powerful framework to bridge microscopic pore properties with macroscopic reservoir behavior. Numerous studies have demonstrated that fractal dimensions correlate with porosity, permeability, and specific surface area, making them indispensable tools for evaluating reservoir quality [10,11,12]. Beyond simple mono-fractal models, multi-fractal statistics and wavelet-based methods allow decomposition of complex pore structures across multiple scales [13,14]. These approaches have extended fractal analysis from static rock characterizations to broader geoscientific contexts, encompassing unconventional reservoir systems, soil structures, and engineered porous materials [15,16,17].

Despite these advances, most existing studies remain limited to static representations of fractal dimensions, neglecting the dynamic evolution of pore structures under fluid-rock interactions and reservoir stress conditions. This limitation obscures the spatiotemporal heterogeneity of pore networks, particularly in systems undergoing continuous diagenesis, hydrocarbon production, or CO₂ injection. Static fractal descriptors fail to capture time-dependent scaling laws or the feedback between pore-scale processes and reservoir-scale behavior [18,19,20]. Consequently, the lack of dynamic fractal frameworks hampers predictive modeling of multi-phase flow, storage capacity, and long-term reservoir performance.

In this work, we present a Dynamic Fractal Dimension Modeling (DFDM) framework that integrates image analysis, wavelet-based fractal methods, and structural fractal geometry to quantify the spatiotemporal evolution of pore networks. Unlike conventional static approaches, our method reveals dynamic scaling laws and demonstrates how fractal dimensions evolve under reservoir-relevant conditions. This framework advances the theoretical understanding of fractal pore dynamics while providing a predictive tool with direct implications for unconventional hydrocarbon recovery, geological CO₂ sequestration, and multi-phase flow in porous media.

2. Modeling

We propose a Dynamic Fractal Dimension Modeling (DFDM) framework that explicitly incorporates temporal variations into fractal descriptors. The central idea is that pore networks evolve as a combined effect of intrinsic reservoir stress-sensitivity and extrinsic fluid-rock interactions. By doing so, DFDM establishes a dynamic, process-oriented bridge between pore-scale complexity and reservoir-scale properties, thereby laying the foundation for the modeling framework presented in this section.

2.1. Fundamentals of Fractal Geometry and Scaling Laws

Fractal geometry provides a powerful mathematical framework for characterizing the irregularity and heterogeneity of pore networks across scales. Unlike Euclidean measures, the fractal dimension D quantifies scale-dependent complexity. The classical box-counting dimension D is defined as:

$$D=\underset{ϵ\to 0}{\mathit{lim}}{\displaystyle \frac{\mathit{log}N\left(ϵ\right)}{\mathit{log}\left(1/ϵ\right)}}$$

(1)

Here, $N\left(ϵ\right)$ denotes the number of covering boxes of size length $\mathsf{ϵ}$ required to cover the pore space. This approach has been widely implemented in tools such as ImageJ with the FracLac plugin, offering a fast and intuitive way to characterize pore morphology. However, such static approaches provide only a single-valued fractal dimension and cannot capture the temporal evolution of pore structures. In reservoir systems, pore networks evolve continuously under stress perturbations, fluid-rock interactions, and diagenetic processes, requiring a more advanced approach.

To address these limitations, wavelet-based fractal analysis has been introduced (see Section 2.3). Unlike the box-counting method, wavelets possess strong time-frequency localization, making them particularly effective in analyzing non-stationary signals, local heterogeneity, and dynamic evolution. This establishes the theoretical basis for the development of the Dynamic Fractal Dimension Model (DFDM) presented later.

2.2. Probability Distribution of Pore Networks

Beyond scaling descriptors, the stochastic nature of pore geometry can be quantified via probability distribution functions (PDFs). The Weibull distribution provides a flexible representation of pore radius statistics:

$$p\left(r\right)={\displaystyle \frac{m}{r_0}}{\left({\displaystyle \frac{r}{r_0}}\right)}^{m-1}exp\left[-{\left({\displaystyle \frac{r}{r_0}}\right)}^m\right]$$

(2)

Here,$m$ is the Weibull modulus reflecting heterogeneity, and $r_0$is the scale parameter representing characteristic pore size. For $m\to \infty$, the distribution converges to a delta function corresponding to homogeneous pores;conversely, small $m$values indicate wide pore size dispersion. This statistical representation provides a baseline for linking probability models with fractal scaling laws. The cumulative distribution function (CDF) is expressed as:

$$F\left(r\right)=1-exp\left[-{\left({\displaystyle \frac{r}{r_0}}\right)}^m\right]$$

(3)

Which can be directly linked to porosity $\varphi$ by integrating over pore volumes:

$$\varphi ={\int }_0^{\infty }{\displaystyle \frac{4}{3}}\pi r^{3}p\left(r\right)dr$$

(4)

To reconcile statistical and fractal descriptions, we impose the condition that the pore number at scale $r$ follows a power law:

$$N\left(r\right)\propto r^{-D}$$

(5)

By comparing with the Weibull PDF, the effective fractal dimension can be expressed as:

$$D_{eff}=-{{\displaystyle \frac{d\mathit{log}p\left(r\right)}{\partial \mathit{log}r}}=m\left({\displaystyle \frac{r}{r_0}}\right)}^m$$

(6)

Thus, the effective fractal dimension${\mathrm{D}}_{\mathrm{e}\mathrm{f}\mathrm{f}}$is not constant but evolves dynamically as pore size distribution parameters $\left(m,r_0\right)$change under stress and fluid interactions. This formulation explicitly reflects scale-dependent irregularity and bridges probabilistic and fractal models.

Limiting Cases

Homogeneous Reservoir ($m\to \infty$):$D_{eff}\to \infty$, indicating near-Euclidean uniform pore geometry.

Highly Heterogeneous Reservoir (($m\to 1$):$D_{eff}$exhibits strong dependence on pore size $r$, reflecting multi-scale irregularity.

Dynamic Evolution: If $r_0{=r}_0\left(t\right)$evolves with compaction, $D_{eff}\left(t\right)$explicitly tracks spatiotemporal heterogeneity.

2.3. Dynamic Fractal Dimension Model (DFDM)

The proposed Dynamic Fractal Dimension Modeling (DFDM) framework establishes a systematic workflow to bridge experimental pore-scale characterization with reservoir-scale predictive modeling. Unlike static fractal approaches (e.g., Image J/FracLac), DFDM explicitly incorporates temporal variability and dynamic interactions, making it particularly suited for stress-sensitive tight reservoirs.

The framework integrates five interconnected modules (Figure 1): (i) acquisition of pore-scale digital images, (ii) multi-scale decomposition using wavelet transforms, (iii) derivation of dynamic fractal dimensions through structural geometry and time-dependent corrections, and (iv) validation against experimental measurements. This hierarchical workflow transforms raw experimental data into dynamic descriptors of pore evolution, ultimately enabling numerical simulations of multi-phase flow and storage processes.

2.3.1. Governing Principles and Dynamic Scaling Law

For pore structures imaged at resolution $ϵ$, this definition of fractal dimension naturally extends to describe the scale-occupation relation over time. Within DFDM, $D$ is allowed to vary dynamically with the physical state of the reservoir, reflecting the coupled influence of porosity change $∆\varphi$,stress perturbation $∆\sigma$and fluid transport volume $∆V$:

$$D\left(t+∆t\right)=D\left(t+∆t\right)+\mathcal{F}\left(∆\varphi ,∆\sigma ,∆V\right)$$

(7)

This equation expresses the change in the fractal dimension$D\left(t\right)$over a small time interval $∆t$,where $\mathcal{F}\left(∆\varphi ,∆\sigma ,∆V\right)$characterizes the nonlinear coupling, establishing a process-based link between pore network evolution and fractal descriptors.

Formally, the dynamic fractal dimension is defined from the time-dependent occupation number $N\left(ϵ,t\right)$ as the local scaling exponent:

$$D\left(t\right)=\underset{ϵ\to 0}{\mathit{lim}}{\displaystyle \frac{\partial \mathit{log}N\left(ϵ,t\right)}{\partial \mathit{log}\left(1/ϵ\right)}}$$

(8)

This expression captures the evolution of pore structures at different scales $ϵ$ over time $t$.To account for stochastic micro-structural variability and non-stationary topology, we generalize the classical power law to:

$$N\left(ϵ,t\right)=C\left(t\right){ϵ}^{-D\left(t\right)}+\eta \left(ϵ,t\right)$$

(9)

where $C\left(t\right)$is a normalization associated with porosity and$\eta \left(ϵ,t\right)$accounts for random micro-structural variability, This decomposition cleanly separates deterministic evolution $\left(capturedbyD\left(t\right)\right)$ from stochastic variability.

Wavelet analysis provides the necessary multi-resolution machinery to track these dynamics. The continuous wavelet transform $W_f\left(a,b\right)$ $\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{e}-\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{s}$ $a$ and positions $b$.

$$W_f\left(a,b\right)={\int }_{-\infty }^{\infty }f\left(x\right){\psi }^{*}\left({\displaystyle \frac{x-b}{a}}\right)dx$$

(10)

And its energy spectrum$E\left(a\right)$ follows a fractal scaling law with Hurst exponent $H$：

$$E\left(a\right)\propto a^{2H+1}$$

(11)

The Hurst exponent $H$ quantifies the long-term memory of pore structure evolution, reflecting the self-similar behavior of the pore network across different scales.

In 2D slices (or 3D volumes with the appropriate embedding dimension $E$), the relation $D\left(t\right)=E-H\left(t\right)$ links spectral roughness to geometric complexity: larger$H\left(t\right)$ implies smoother, more persistent connectivity (smaller $D$); smaller $H\left(t\right)$ implies irregular, fracture-dominated patterns (larger $D$).

Under static equilibrium, $N\left(ϵ,t\right)$is time-invariant and$\partial N\left(ϵ,t\right)=0$, then the fractal dimension $D\left(t\right)$ reduces to a constant value $D_0$. However, under dynamic conditions (such as during compaction, fracture development, or fluid injection), the fractal dimension $D$evolves as a function of both time and pore structure changes, reflecting the dynamic coupling between mechanical compaction, fluid flow, and diagenesis. The evolution of pore networks in tight reservoirs arises from two dominant processes: mechanical compaction and stress redistribution, reducing pore volume; Fluid-rock interactions, dissolving pores or fractures at multiple scales. We therefore define the dynamic fractal dimension$D\left(t,x\right)$ as a spatiotemporal field variable, dependent on both time $t$ and spatial coordinate $x$:

$$D\left(t,x\right)=\underset{ϵ\to 0}{\mathit{lim}}{\displaystyle \frac{\partial \mathit{log}N\left(ϵ,t,x\right)}{\partial \mathit{log}\left(1/ϵ\right)}}$$

(12)

Here, $N\left(ϵ,t,x\right)$ denotes the local occupation number of pores at resolution $ϵ$, time $t$, and position $x$. Building on this spatiotemporal definition, the classical scaling laws must be generalized to incorporate temporal variability, motivating a Dynamic Scaling Law in which time-dependent coefficients capture the evolving effects of stress, fluid flow, and diagenesis.

$$N\left(ϵ,t\right)=C_0\left(t\right)e^{-at}{ϵ}^{-D\left(t\right)}+\eta \left(ϵ,t\right)$$

(13)

Here, $C_0$represents the initial scaling constant linked to porosity at$t=0$,$a$ is the compaction rate parameter that governs pore size evolution. By differentiating the generalized power-law scaling equation with respect to time, we derive the dynamic scaling equation, which links the time-dependent evolution of the fractal dimension to the temporal variation in the pore network structure:

$${\displaystyle \frac{\partial D}{\partial t}}=-{\displaystyle \frac{1}{\mathit{log}\left(1/ϵ\right)}}{\displaystyle \frac{{\partial }^2\mathit{log}N\left(ϵ,t\right)}{\partial t\partial \mathit{log}ϵ}}$$

(14)

This equation quantifies how the internal reorganization of the pore network modulates the rate of change of fractal dimension.

2.3.2. Wavelet Spectral Estimation

The wavelet energy spectrum of pore images evolves with time as:

$$E\left(a,t\right)=\int {\left|W_f\left(a,b,t\right)\right|}^{2}db$$

(15)

This equation reflects the redistribution of pore-scale heterogeneity across scales, allowing dynamic tracking of structural evolution, where a is scale, b is position, and t is time. Empirically,

$E\left(a,t\right)\propto a^{2H\left(t\right)+1},$ enabling a direct time-resolved estimate of the Hurst exponent and, via$D\left(t\right)=E-H\left(t\right),$ of the dynamic fractal dimension. Physically, higher $H\left(t\right)$indicates more persistent, better-connected structures; lower$H\left(t\right)$indicates increased irregularity or the emergence of fractures and channels.

To incorporate stress and flow effects in a tractable way for reservoir-scale simulations, we approximate the time-varying fractal dimension$D\left(x,t\right)$inferred from wavelet spectra by a reduced-order phenomenological expression:

$$D\left(x,t\right)=D_0-\alpha {\displaystyle \frac{t}{t_{max}}}-\beta {\displaystyle \frac{x}{L}}+\gamma \sin \left(\pi {\displaystyle \frac{x}{L}}\right)e^{-\delta t/t_{max}}$$

(16)

Here, $D_0$ denotes the static baseline fractal dimension (≈2.57 in this study), $\alpha$ and $\beta$ quantify temporal and spatial attenuation, respectively, $\gamma$represents the local oscillations, and$\delta$damps those oscillations. This surrogate preserves the dominant multi-scale trends revealed by wavelets while keeping the subsequent numerical model computationally efficient.

2.3.3. Coupled Fluid-Solid Interaction

Beyond spectral estimation, pore structure dynamics are also strongly influenced by external driving forces, particularly the coupling between stress perturbations and fluid invasion. We define an effective fractal dimension $D_{eff}\left(t\right)$ governed by the competition between mechanical compaction (stress tensor ${\sigma }_{ij}$) and flow-driven enhancement (pressure gradient $\nabla p$):

$${\displaystyle \frac{d{D}_{eff}}{dt}}=-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖$$

(17)

Here,$\beta$quantifies stress sensitivity (volumetric stress promotes pore closure and structural simplification), while$\gamma$ quantifies fluid-induced dissolution/scouring and the opening of new flow channels.

Diagnostic regimes:

Pure mechanical compaction ($‖\nabla p‖=0$):${\displaystyle \frac{d{D}_{eff}}{dt}}=-\beta tr\left({\sigma }_{ij}\right)<0$,fractal dimension decreases monotonically as pores collapse under stress.

Pure fluid dissolution $\left({\sigma }_{ij}=0\right):{\displaystyle \frac{d{D}_{eff}}{dt}}=\gamma ‖\nabla p‖>0$,fractal dimension increases due to flow-driven pore opening and dissolution.

Steady-state equilibrium $\beta tr\left({\sigma }_{ij}\right)=\gamma ‖\nabla p‖:{\displaystyle \frac{d{D}_{eff}}{dt}}=0$,pore system reaches a dynamic equilibrium where compaction and dissolution exactly balance, maintaining a stable fractal state.

This formulation indicates that compaction processes generally decrease fractal complexity, whereas fluid-solid interactions driven by strong pressure gradients enhance it, highlighting the competing roles of stress and flow in governing pore evolution, providing a first-order closure linking fractal descriptors to measurable geomechanical and hydraulic parameters.

2.3.4. Dynamic Porosity-Permeability Relationship

To connect fractal dynamics to transport properties, we express porosity$\varphi \left(t\right)$and permeability $k\left(t\right)$as functions of$D\left(t\right).$Following Katz-Thompson-type scaling, permeability obeys a fractal-dependent power law:

$${\displaystyle \frac{d{D}_{eff}}{dt}}=-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖$$

(18)

The equation describes how permeability evolves over time as a function of porosity, with the fractal dimension $D\left(t\right)$dictating the scaling relationship between these two properties. This is based on the assumption that changes in the pore network structure, as characterized by$D\left(t\right)$, significantly influence the connectivity of the pore spaces, which in turn impacts the flow characteristics of the reservoir.

Mechanistically, decreasing$D\left(t\right)$under compaction generally correlates with declining permeability, while increasing$D\left(t\right)$under fluid-induced reorganization correlates with enhanced permeability through improved connectivity. Tracking the evolution of $D\left(t\right)$therefore provides a predictive handle on how reservoir processes (compaction, dissolution, fracture activation) reshape the permeability-porosity relation over time.

2.4. Governing Equation for Spatiotemporal Fractal Evolution

2.4.1. Field Representation

We consider the fractal dimension as a continuous scalar field $D\left(x,t\right)$ defined over the reservoir domain $\mathsf{\Omega }$. Its dynamics are driven by competing processes: Compaction-driven reduction, caused by mechanical stress and pore closure; Fluid-driven enhancement, arising from dissolution, fracturing, and reorganization; Diffusive smoothing, redistributing pore-scale heterogeneity across scales. Thus, the general governing PDE takes the form:

$${\displaystyle \frac{\partial D}{\partial t}}=\nabla \bullet \left(\kappa \nabla D\right)-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖+\xi \left(x,t\right)$$

(19)

Here, $\kappa$is the fractal diffusivity (smoothing coefficient), the diffusive term $\nabla \bullet \left(\kappa \nabla D\right)$smooths local fluctuations in pore-scale heterogeneity; $\beta$is the stress-compaction coefficient, the compaction term$-\beta tr\left({\sigma }_{ij}\right)$accounts for stress-driven pore closure; $\gamma$is the fluid enhancement coefficient, the fluid-driven enhancement term $\gamma ‖\nabla p‖$ reflects dissolution and flow-induced reorganization of pore networks. The stochastic forcing term $\xi \left(x,t\right)$captures random fracture nucleation, mineral dissolution, and other unresolved micro-scale processes.

This PDE naturally embeds the dynamic porosity-permeability relationship (Section 2.3.4) within a spatiotemporal framework. Specifically, the solution of $D\left(x,t\right)$ governs the evolution of porosity $\varphi \left(x,t\right)$ and permeability $k\left(x,t\right)$ across the reservoir domain, thereby coupling micro-scale fractal dynamics with macro-scale transport properties.

Porosity $\varphi \left(x,t\right)$ evolves in tandem with fractal dimension, a nonlinear coupling is assumed:

$${\displaystyle \frac{\partial \varphi }{\partial t}}=-\alpha \varphi +\lambda f\left(D\left(x,t\right)\right)$$

(20)

Here, $\alpha$ is a compaction rate and $f\left(D\right)$is a monotonic function, e.g.,$f\left(D\right){=\varphi }_{max}{\left({\displaystyle \frac{D}{D_c}}\right)}^m$.

Similarly,Permeability $k\left(x,t\right)$ follows a fractal-dependent scaling (cf. Katz-Thompson):

$$k\left(x,t\right)\propto {\varphi \left(x,t\right)}^{{\displaystyle \frac{D\left(x,t\right)}{3-D\left(x,t\right)}}}$$

(21)

This formulation highlights the bidirectional feedback between fractal dynamics and reservoir properties: fractal evolution modifies porosity and permeability, while in turn stress and fluid-rock interactions drive changes in fractal geometry.

2.4.2. Limiting Cases of the PDE

The general PDE admits several limiting forms under specific physical conditions, which provide physical intuition for the role of diffusion, stress, and flow.

Case 1: Pure Diffusion $\left(\beta =\gamma =0\right),{\displaystyle \frac{\partial D}{\partial t}}=\kappa {\nabla }^{2}D$.

This case corresponds to the redistribution of pore-scale heterogeneity without net gain or loss of pore volume, analogous to a heat equation-pore complexity homogenizes over time.

Case 2: Stress-Dominated $\left(\gamma =0\right),{\displaystyle \frac{\partial D}{\partial t}}=\kappa {\nabla }^{2}D-\beta tr\left({\sigma }_{ij}\right).$Under this regime, compaction continuously decreases pore connectivity, fractal dimension decreases steadily with compaction monotonically.

Case 3: Flow-Dominate$\left(\beta =0\right),{\displaystyle \frac{\partial D}{\partial t}}=\kappa {\nabla }^{2}D+\gamma ‖\nabla p‖.$Here, Dissolution and flow-driven reorganization promote the opening of new pathways, leading to an increase in fractal complexity.

These limiting cases highlight the fundamental physical regimes captured by the governing PDE, ranging from pure diffusion to stress- and flow-dominated evolution.

2.5. Numerical Implementation Strategy

2.5.1. Discretization

We discretize the governing PDE using a finite element/finite volume method, depending on the geometry of reservoir samples. The weak formulation is:

$${{\int }_{\Omega }\omega {\displaystyle \frac{\partial D}{\partial t}}d\Omega +\int }_{\Omega }\kappa \nabla \omega \bullet \nabla Dd\Omega {=\int }_{\Omega }\omega \left(-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖+\xi \right)d\Omega$$

(22)

Here, $\mathsf{\omega }$ is a test function. In practice, dynamic fractal dimension calculations can be implemented using both box-counting and wavelet-based methods. MATLAB routines and ImageJ plugins (e.g., FracLac) enable efficient fractal estimation and numerical stability. This combined strategy highlights the advantage of wavelet transforms in resolving multi-scale dynamics, whereas box-counting offers a robust baseline for static scaling analysis.

2.5.2. Multi-scale Up-scaling and Stability Constraints

Since reservoir processes span nano-meters to meters, we adopt a two-scale strategy to ensure both physical and numerical stability.

Micro-scale (cast thin-section images): Fractal evolution is extracted at representative elementary volumes (REV), typically using high-resolution pore-scale imaging.

Macro-scale (reservoir grid): Effective fields of $D\left(x,t\right)$,$\varphi \left(x,t\right)$and$k\left(x,t\right)$are obtained by homogenization and transferred into reservoir-scale simulations.

To guarantee well-posed numerical evolution during up-scaling, additional stability constraints are imposed:

Positivity Constraint: $0<D\left(x,t\right)<3$.

Energy Dissipation Constraint: ${\displaystyle \frac{d}{dt}}{\int }_{\Omega }{D}^{2}d\Omega \le 0$.

These constraints prevent nonphysical divergence and guarantee monotonic energy decay, ensuring robust numerical solutions during up-scaling.

2.6. Validation and Experimental Bench-Marking

2.6.1. Experimental Data and Calibration

Validation of the DFDM framework relies on high-resolution experimental datasets:

Case Thin-Section Images: Provide high-resolution 2D pore maps under stress and fluid perturbations, enabling direct assessment of spatiotemporal fractal evolution.

Displacement experiments: Capture multi-phase flow behavior and fluid-rock interactions, serving as benchmarks for evaluating the predictive capability.

Calibration is performed by computing spatial correlation functions:

$$\mathrm{C}\left(r,t\right)=〈\delta D\left(x,t\right)\delta D\left(x+r,t\right)〉,\delta D=D-\overline{D}$$

(23)

Ensuring that the simulated pore structures and flow dynamics reproduced the experimental observations from these datasets. Importantly, cast thin-section imaging constrains the structural heterogeneity of the pore system, while displacement experiments constrain the dynamic flow behavior, making them complementary to the DFDM framework. This integration enhances the reliability and robustness of the subsequent simulations.

2.6.2. Model Validation and Bench-marking

Validation proceeds through comparison between DFDM predictions and experimental measurements:

Static Validation: At $\mathrm{t}=0$, the predicted baseline fractal dimension $D\left(0\right)$was benchmarked against values derived from case thin-section images (Figure 2). The close agreement (deviation <5%) ensured correct initialization of DFDM simulations.

Dynamic Validation: The temporal evolution of $D\left(t\right)$predicted by DFDM was compared against experimental curves obtained from CO₂ flooding through long cores under controlled injection conditions. The error metric was quantified as:

$${\epsilon }_D\left(t\right)={\displaystyle \frac{\left|D_{sim}\left(t\right)-D_{exp}\left(t\right)\right|}{D_{exp}\left(t\right)}}$$

(24)

With acceptable agreement defined as ${\epsilon }_D\left(t\right)<10\%$ across multiple time scales.

Although explicit boundary-condition-limiting tests (e.g., pure compaction, pure dissolution) were not performed, the CO₂ flooding experiments inherently encapsulated competing effects of stress-induced compaction and fluid-driven dissolution. This indirect validation demonstrated the model’s ability to capture spatiotemporal evolution of pore networks under realistic reservoir conditions.

Figure 2 illustrates the image-to-fractal workflow, where thin-section micro-graphs were processed into binary pore-matrix representations and analyzed via box-dimension regression, yielding baseline fractal descriptors for DFDM initialization and subsequent validation.

In this context, image-based analysis provides an intuitive bridge between experimental observation and model validation. As illustrated in Figure 2, the workflow begins with thin-section micro-graphs of the rock matrix (a), which are subsequently converted into binary pore-matrix representations (b). These segmented images are then subjected to box-dimension regression using the MATLAB interface (c), yielding quantitative fractal dimensions. This sequence explicitly demonstrates how pore-scale heterogeneity, captured through imaging, is rigorously translated into mathematical descriptors that can be bench-marked against the predictions of the Dynamic Fractal Dimension Model (DFDM).

2.6.3. Statistical Confidence

The robustness of DFDM was further established by assessing statistical confidence across experimental datasets. Quantitative agreement between model predictions and experimental measurements was evaluated using the coefficient of determination:

$$R^2=1-{\displaystyle \frac{{\sum }_t{\left(D_{sim}\left(t\right)-D_{exp}\left(t\right)\right)}^2}{{\sum }_t{\left(D_{exp}\left(t\right)-{\overline{D}}_{exp}\right)}^2}}$$

(25)

A predictive model requires $R^2>0.9$. The obtained values exceeded this threshold, confirming that DFDM reliably reproduces both structural constraints from thin-section analysis and flow constraints from CO₂ flooding experiments. By explicitly integrating temporal variability, fluid-solid coupling, and spatiotemporal PDEs, DFDM provides a rigorous predictive framework for pore-structure evolution in stress-sensitive reservoirs, with direct implications for CO₂ storage and enhanced hydrocarbon recovery.

3. Results and Discussion

The proposed Dynamic Fractal Dimension Modeling (DFDM) framework provides a systematic approach to quantify the spatiotemporal evolution of pore structures in tight reservoirs, bridging pore-scale imaging with predictive reservoir-scale simulations. By explicitly incorporating temporal variability into fractal characterization, DFDM represents a paradigm shift from static pore-scale descriptors to dynamic, predictive modeling of reservoir evolution.

In this section, we present the numerical predictions, experimental validations, and mechanistic insights derived from DFDM, followed by a critical discussion of model robustness, limitations, and implications for reservoir engineering. To assess the predictive capacity of the framework, model outputs were directly compared with laboratory measurements, sensitivity analysis with alternative parameter sets underscores the model’s adaptability to varying geological conditions.

Together, these results confirm that DFDM not only advances the theoretical understanding of fractal pore dynamics but also provides a practical tool for predicting reservoir behavior in unconventional hydrocarbon recovery and CO₂ sequestration. The simulation results form the basis for the following discussion, including spatiotemporal evolution, mechanistic insights, and model validation and limitations.

3.1. Spatiotemporal Evolution of Fractal Dimension

The spatiotemporal evolution of fractal dimensions under CO₂ flooding reflects the coupled effects of pore-scale heterogeneity, stress sensitivity, and multi-phase flow. To capture this process, a stepwise framework is employed, beginning with static baseline derived from thin-section imaging, followed by experimental observation of pressure propagation, and ultimately dynamic prediction through the DFDM approach. Figure 3, Figure 4 and Figure 5 collectively illustrate this progression: the static fractal baseline (Figure 3), the experimentally observed displacement stages (Figure 4), and the predicted dynamic evolution of fractal dimension (Figure 5).

The static fractal dimension obtained from thin-section images provides a necessary reference for evaluating the dynamic evolution of pore structures under CO₂ flooding. As illustrated in Figure 3, the box-counting regression applied to binary micro-graphs yields a fractal dimension of approximately $D_0\approx 2.57$.

To further investigate displacement behavior, pressure-propagation experiments were carried out along a long-core sample during CO₂ injection. The measured pressure responses at nine monitoring points reveal systematic changes over time. As shown in Figure 4, three distinct stages can be identified: (i) the mixed-phase formation stage, (ii) the mixed-phase displacement stage, and (iii) the gas post-breakthrough stage. These transitions mark the progressive redistribution of pore pressure and reflect the underlying stress-sensitive heterogeneity of the rock matrix.

To further interpret these experimental findings, DFDM was applied to predict the spatiotemporal distribution of fractal dimension $D\left(x,t\right)$,the results are presented in Figure 5. The contour map demonstrates a systematic evolution from higher values near the injection point toward lower values after breakthrough as CO₂ flooding proceeds. Notably, the baseline dimension $D_0\approx 2.57$ serves as a transition threshold: values above this baseline correspond to relatively well-connected pore networks during early displacement, while reductions below$D_0$mark the onset of channelization and non-uniform flow. The predicted dynamic range (1.5-2.8) aligns with experimental observations and remains physically consistent with pore-scale heterogeneity and stress-induced permeability alterations.

3.2. Model Validation and Predictive Robustness

The predictive capability of the Dynamic Fractal Dimension Model (DFDM) must be validated to ensure its applicability under reservoir-relevant conditions. While the static fractal dimension ($D_0\approx 2.57$) has been experimentally determined and serves as the baseline descriptor of pore structures, direct measurements of dynamic fractal dimensions are not available. To address this limitation, the validation strategy adopted here relies on three complementary aspects: (i) consistency of predicted trends with experimental pressure responses and pore-network evolution, (ii) alignment of the initial condition with the experimentally derived static fractal dimension, and (iii) agreement with scaling behaviors reported in previous studies.

As illustrated in Figure 6, the DFDM framework predicts a progressive decline in the fractal dimension $D\left(V\right)$ with increasing cumulative CO₂ injection volume. Two parameter sets are shown to represent the range of model variability (solid and dashed lines). The experimental static fractal dimension ($D_0\approx 2.57$) is included as a reference point at the onset of injection, serving as the critical anchor for validation. Although direct dynamic measurements are absent, the predicted reduction of D is physically consistent with pore-network simplification, compaction-induced stress redistribution, and non-uniform flow patterns during CO₂ flooding. This agreement underscores the physical plausibility of DFDM predictions.

Beyond validation, A systematic sensitivity analysis was performed to evaluate the robustness of DFDM predictions under parameter perturbations. Variations in the scaling exponent, stress-coupling coefficient, and initial pore heterogeneity were introduced to examine their influence on the predicted fractal dimension $D\left(t\right)$. Results (Figure 7) demonstrate that while parameter changes shift the detailed trajectory of$D\left(t\right)$, the overall predictive fidelity and mechanistic stability of the framework are preserved.

3.3. Mechanistic Insights: Stress Sensitivity and Flow-Driven Processes

The mechanistic analysis reveals the fundamental competition between stress-induced compaction (reducing D and permeability) and flow-driven dissolution/enhancement (increasing D and improving connectivity), As illustrated in Figure 8, these dual regimes delineate the operational bounds of reservoir evolution:

Compaction-dominated regime: Pore contraction reduces connectivity, leading to permeability loss.

Flow-driven regime: Fluid injection enlarges pore structures, accelerating flow pathways.

By explicitly resolving these dual processes, DFDM provides a mechanistic closure between fractal descriptors and stress-flow interactions, extending beyond descriptive statistics toward predictive physical modeling of pore-network evolution.

3.4. Comparative Performance and Model Limitations

While DFDM exhibits substantial predictive capability, several limitations warrant critical evaluation. First, the computational demand remains high, limiting large-scale applications without advanced solvers. Second, the predictive accuracy of DFDM is strongly dependent on the representativeness of input data, reconstructions derived from 2D thin-section images, for instance, may under-sample sub-resolution heterogeneities that are critical for transport dynamics, thereby introducing scale-dependent biases. Third, the model’s sensitivity to stress, fracture-pore coupling coefficients, and boundary conditions introduces non-negligible uncertainties when extrapolating beyond controlled laboratory settings.

Nonetheless, the hierarchical architecture of DFDM renders it highly extensible. Integration with machine-learning-assisted parameter estimation, physics-informed neural networks, or hybrid multi-scale solvers could substantially reduce computational demands and enhance generalizability. Despite its limitations, DFDM offers a robust, physically grounded pathway for bridging pore-scale heterogeneity with reservoir-scale performance.

4. Conclusions

This study introduces a Dynamic Fractal Dimension Modeling (DFDM) framework to quantify the spatiotemporal evolution of pore structures in tight reservoirs. By integrating high-resolution thin-section imaging, CO₂ long-core displacement experiments, wavelet-based fractal methods, and structural fractal geometry, DFDM advances beyond conventional static fractal approaches to establish a dynamic, time-resolved representation of pore heterogeneity under reservoir-relevant conditions.

The results show that DFDM successfully captures critical transitions in pore complexity during CO₂ flooding, linking stress-induced compaction and flow-driven reorganization with measurable changes in fractal dimension. Model validation against laboratory experiments demonstrates strong agreement with pressure responses, while sensitivity analysis confirms robustness across varying geological settings.

The framework provides direct implications for reservoir engineering: (i) predicting permeability loss and optimizing injection strategies in enhanced oil recovery, (ii) evaluating capacity and long-term stability in geological CO₂ storage, and (iii) extending to geothermal systems, aquifers, and engineered porous media where heterogeneity governs multi-phase flow.

In summary, DFDM introduces a dynamic and predictive approach to pore-scale heterogeneity, bridging imaging, fractal geometry, and stress-coupled transport. This methodological innovation not only advances theoretical understanding of fractal pore dynamics but also establishes a practical tool for energy recovery, carbon sequestration, and sustainable subsurface resource management. DFDM unifies pore-scale imaging and reservoir-scale modeling, offering a robust framework to quantify and predict dynamic heterogeneity in subsurface energy systems.