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 CO₂ storage efficiency. However, the intrinsic heterogeneity and complexity of these formations continue to pose major challenges for accurate reservoir characterization [1,2]. Shales, coals, and carbonate rocks exhibit multi-scale pore systems ranging from nanometre-scale micro-pores to interconnected micro-fractures, forming a hierarchical structure with significant spatial variability [3,4]. This heterogeneity governs adsorption, permeability, and capillary behavior, thereby influencing reservoir quality and production performance [5,6].

Recent advances in experimental and numerical techniques have greatly improved our understanding of pore systems. Methods such as low-pressure gas adsorption, NMR relaxometry, and X-ray micro-CT imaging have revealed complex pore geometries and connectivity patterns across scales [7,8]. Quantitative characterization of these pore networks has demonstrated that their morphology and surface roughness often obey fractal scaling laws, implying that fractal geometry provides a powerful mathematical framework for describing pore heterogeneity [9,10,11]. By linking geometric irregularity to measurable physical properties, fractal analysis enables the quantification of porosity, specific surface area, tortuosity, and permeability in unconventional reservoirs [12,13,14]. Building on this foundation, recent studies have refined fractal and multi-fractal descriptors to better capture heterogeneity across scales [15,16], enhancing our ability to relate micro-structural features to macroscopic transport behavior.

The concept of fractal dimensions has thus become an indispensable tool in reservoir petrophysics. Numerous studies have confirmed that fractal dimensions correlate strongly with reservoir parameters such as porosity, pore-size distribution, and permeability [17,18,19]. Notably, recent studies have demonstrated the applicability of pore-type-dependent fractal analysis to shale reservoirs [20] and advanced multi-fractal scaling for dynamic heterogeneity characterization [21,22]. Furthermore, multi-fractal analysis and wavelet-based methods have been successfully applied to distinguish the contributions of different pore scales and to describe the dynamic evolution of pore morphology [23,24,25]. Such approaches provide insights into the scaling and self-similarity of pore systems, revealing the hierarchical nature of storage and flow mechanisms within tight reservoirs [26,27].

Despite these advances, traditional fractal approaches typically treat pore geometry as static, neglecting the time-dependent evolution of pore structures under stress, temperature, and fluid-rock interaction [28,29]. This limitation hampers predictive modeling because actual reservoirs experience continuous structural evolution due to compaction, dissolution-precipitation reactions, and multi-phase flow [30,31]. Dynamic fractal models have been proposed to capture this temporal evolution [32], thereby establishing a conceptual bridge between micro-structural evolution and macroscopic transport behavior [33,34,35]. Such models enable the interpretation of hysteresis and permeability changes by quantifying how the fractal dimension evolves with stress and saturation [36].

Several researchers have further embedded fractal descriptors into digital rock physics and flow simulations. By coupling micro-CT imaging with lattice Boltzmann or pore-network models, fractal-based parameters can reproduce non-linear flow, capillary hysteresis, and wettability-dependent behavior [37]. Nevertheless, challenges remain in establishing a unified, physics-based dynamic fractal framework that captures the coupling between mechanical deformation, fluid saturation, and evolving pore topology [38,39]. Addressing this gap requires integrating time-lapse observations, scale-localized analysis, and constitutive relations that embed fractal descriptors within transport equations [40,41]. Recent studies have begun to explore this by tracking fractal dimension evolution during stress cycling and chemical alteration, revealing distinct signatures of structural versus topological transformations [42,43]. These insights highlight the potential of fractal geometry to serve as a bridge between micro-structure and macro-scale performance, enabling the formulation of data-efficient, physics-informed models for unconventional reservoir development.

In this context, the present study aims to develop a dynamic fractal characterization framework that couples evolving pore structure with transport and mechanical properties. By quantifying the temporal trajectory of fractal parameters under varying reservoir conditions, this framework provides a comprehensive link between pore-scale morphology and macroscopic flow behavior. Ultimately, it seeks to advance our understanding of how structural heterogeneity, fractal scaling, and dynamic evolution jointly control the performance of tight reservoirs under realistic subsurface environments.

In contrast to previous multi-fractal or wavelet-based approaches that statistically track heterogeneity (e.g., in shale or carbonate systems) and to stress-dependent fractal permeability models that remain static in time, the Dynamic Fractal Dimension Model (DFDM) introduces a physics-informed, time-dependent framework for pore-structure evolution.

The key technical innovations are:

(i) formulation of a partial differential equation for the spatiotemporal evolution of the fractal dimension $D\left(x,t\right)$ that integrates diffusive, stress-induced, and flow-driven mechanisms;

(ii) establishment of explicit coupled relationships between fractal dimension, porosity, and permeability consistent with the Katz-Thompson scaling [44]; and

(iii) incorporation of stochastic forcing to represent micro-scale restructuring under dynamic stress-flow interactions.

This constitutes, to our knowledge, the first attempt to couple a PDE on $D\left(x,t\right)$with porosity-permeability evolution under CO₂-driven or stress-perturbed flow regimes, thereby extending conventional static fractal descriptions toward predictive modeling of reservoir micro-structure dynamics.

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.The DFDM modeling framework is formulated in terms of explicit inputs, outputs, and underlying assumptions. The input data comprise binarized pore-space images obtained from thin-section characterization, pressure-displacement curves along the core length, and external loading parameters describing boundary stress conditions. The primary output of the model is the spatiotemporal field of the dynamic fractal dimension $D\left(x,t\right)$, which captures the evolution of pore geometry under coupled mechanical and flow processes. In addition, auxiliary properties such as porosity$\varphi (x,t)$ and permeability$k(x,t)$ can be inferred from $D\left(x,t\right)$ through established fractal or empirical relationships. These derived quantities serve to bridge the dynamic fractal analysis with experimentally measurable transport parameters.

To ensure clarity and consistency in notation, all variables and parameters used in the modeling framework are summarized in Table A1 of Appendix A.1.

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{\mathsf{\epsilon }\to 0}{\mathit{lim}}{\displaystyle \frac{\mathit{log}N\left(\mathsf{\epsilon }\right)}{\mathit{log}\left(1/\mathsf{\epsilon }\right)}}$$

(2.1)

where $N\left(\epsilon \right)$ denotes the number of covering boxes of size length $\mathsf{\epsilon }$ covering the pore space,with

$\epsilon$ representing the characteristic length scale of the imaging domain (in m). 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. $D$ values derived from the wavelet-based approach were compared with those obtained from the classical box-counting method. Both methods yielded consistent mean values ($\left|\Delta D\right|<0.05$), while $D_w$exhibited smoother temporal trends and improved stability under noisy or non-stationary conditions. This confirms that the wavelet analysis does not alter the underlying fractal scaling but provides enhanced localization and noise resilience. 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.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]$$

(2.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$$

(2.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}$$

(2.5)

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

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

(2.6)

The detailed derivation of Eq.2.6, including the underlying assumptions, is provided in A.2.1 of Appendix A.2. It is worth noting that Eq.2.6 defines a scale-dependent $D_{eff}\left(r\right)$. In the asymptotic regime ${r\ll r}_0$, where the Weibull distribution approaches a power-law form $p\left(r\right)~r^{m-1}$, $D_{eff}$ reduces to $m$. Therefore, $m$ represents the limiting fractal dimension, while $D_{eff}$ varies with both $r$ and the evolving pore-scale parameters $\left(m,r_0\right)$. Thus, the effective fractal dimension $D_{eff}$ 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

Within DFDM, the fractal dimension D evolves dynamically with the physical state of the reservoir, reflecting the coupled influence of porosity change $∆\varphi$, stress perturbation $∆\sigma$, and fluid transport $∆V$:

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

(2.7)

Here, $\mathcal{F}\left(∆\varphi ,∆\sigma ,∆V\right)$ represents a nonlinear coupling function quantifying the rate and direction of fractal evolution under multi-physical perturbations. This dynamic scaling law provides a process-based link between pore-structure evolution and the time-dependent fractal descriptors. Formally, the dynamic fractal dimension is defined from the time-dependent occupation number $N\left(\epsilon ,t\right)$ as the local scaling exponent:,

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

(2.8)

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

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

(2.9)

where $C\left(t\right)$is a normalization associated with porosity and$\eta \left(\epsilon ,t\right)$accounts for random micro-structural variability, This decomposition cleanly separates deterministic evolution (captured by $D\left(t\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)$ decomposes pore-image signals across scales 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$$

(2.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}$$

(2.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 thin sections reconstructed as gray-level scalar fields embedded in ${\mathbb{R}}^3$, the relationship between the Hurst exponent and fractal dimension follows the standard self-affine form$D\left(t\right)=3-H\left(t\right)$,where $H\left(t\right)$ quantifies the spectral roughness obtained from wavelet-based analysis, and the embedding dimension $E=3$ reflects the gray-level surface representation $z=I\left(x,y\right)$.

Higher $H\left(t\right)$indicates smoother, more persistent connectivity (lower$D\left(t\right)$), while smaller $H\left(t\right)$ corresponds to increased irregularity and rougher pore structures (higher$D\left(t\right)$).

This correspondence can be intuitively interpreted in terms of pore morphology. A higher Hurst exponent $H\left(t\right)$ represents local variations in pore radius are gradual and the structure approaches a quasi-continuous texture, thus yielding a lower effective fractal dimension $D\left(t\right)$. Conversely, a lower $H\left(t\right)$characterized by abrupt multi-scale fluctuations in pore boundaries and higher structural roughness, which correspond to a larger $D\left(t\right)$. This morphological contrast, high H/low D versus low H/high D, provides a physical interpretation of the spectral-geometric linkage in heterogeneous reservoirs. This consistent D-H relation establishes the physical linkage between spectral roughness and morphological complexity in the DFDM framework.

Under static equilibrium, $N\left(\epsilon ,t\right)$is time-invariant and$\partial N\left(\epsilon ,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{\epsilon \to 0}{\mathit{lim}}{\displaystyle \frac{\partial \mathit{log}N\left(\epsilon ,t,x\right)}{\partial \mathit{log}\left(1/\epsilon \right)}}$$

(2.12)

Here, $N\left(\epsilon ,t,x\right)$ denotes the local occupation number of pores at resolution $\epsilon$, 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(\epsilon ,t\right)=C_0\left(t\right)e^{-at}{\epsilon }^{-D\left(t\right)}+\eta \left(\epsilon ,t\right)$$

(2.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/\epsilon \right)}}{\displaystyle \frac{{\partial }^2\mathit{log}N\left(\epsilon ,t\right)}{\partial t\partial \mathit{log}\epsilon }}$$

(2.14)

This equation quantifies how the internal reorganization of the pore network modulates the rate of change of fractal dimension. The detailed algorithm for parameter extraction and stability evaluation is provided in A.2.2 of Appendix A.2.

2.3.2. Wavelet Spectral Estimation

The multiscale analysis of pore heterogeneity was conducted using the FracLab 2.2 toolbox, which integrates continuous wavelet transforms and fractal geometry descriptors. The Morlet wavelet was adopted as the mother function for the continuous wavelet transform owing to its superior time-frequency localization and smooth oscillatory kernel, ensuring accurate characterization of pore-boundary irregularities across scales. Prior to scale decomposition, image noise was mitigated through wavelet shrinkage filtering, where high-frequency coefficients were adaptively thresholded to suppress random noise while preserving geometrical edges and connectivity of the pore network. This procedure effectively enhances the signal-to-noise ratio without compromising multiscale roughness, thus improving the stability of the subsequent fractal estimation.

To incorporate the temporal evolution of pore heterogeneity, spatial coordinates were treated as a pseudo-time variable, reflecting the progressive compaction along the sample core. This spatiotemporal mapping enables tracking of the dynamic Hurst exponent $H\left(t\right)$ and the corresponding effective fractal dimension $D\left(t\right)$, linking structural scaling variability with evolving stress and fluid effects.

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$$

(2.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. For gray-level pore surfaces (embedded in ${\mathbb{R}}^3$), the effective fractal dimension is obtained using $D=3-H$.This formulation is consistent with the self-affine scaling theory of natural surfaces and provides a physically justified link between wavelet-based roughness and fractal geometry.

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 \mathit{sin}\left(\pi {\displaystyle \frac{x}{L}}\right)e^{-\delta t/t_{max}}$$

(2.16)

Here, $D_0$ denotes the static baseline fractal dimension (≈2.57 in this study).

The coefficients $\alpha$ and $\beta$ quantify the temporal and spatial attenuation of the fractal dimension,, respectively (both dimensionless). The term $\gamma \mathit{sin}\left(\pi {\displaystyle \frac{x}{L}}\right)e^{-\delta t/t_{max}}$represents represents small-scale oscillations in local pore geometry, where $\gamma$(dimensionless) controls the oscillation amplitude and local oscillations, and $\delta$(dimensionless) governs its exponential damping over time. This oscillatory term captures the periodic heterogeneity of pore throats and grain contacts observed in thin-section images, while the decay factor$e^{-\delta t/t_{max}}$ensures that these fluctuations gradually vanish as the system approaches a quasi-steady fractal state. Such a reduced-order formulation maintains physical interpretability, while keeping the model computationally efficient for core-scale simulations and readily extensible to reservoir-scale applications.

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‖$$

(2.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:

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

(2.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 $\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 ·\left(\kappa \nabla D\right)-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖+\xi \left(x,t\right)$$

(2.19)

Here, $\kappa$ is the fractal diffusivity (smoothing coefficient), the diffusive term $\nabla ·\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. The numerical discretization and solver configuration for Eq 2.12 detailed in Appendix A.4.

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.

To describe the dependence of porosity $\varphi \left(x,t\right)$on the evolving fractal morphology,porosity evolves in tandem with fractal dimension, a nonlinear empirical coupling relationship proposed [45,46]:

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

(2.20)

where $\alpha$ and$\lambda$denote the stress-induced compaction and structure-driven dilation rates, respectively. The function $f_{\varphi }\left(D\right)$characterizes the influence of fractal geometry on pore connectivity and is expressed as :$f_{\varphi }\left(D\right){=\varphi }_{max}{\left({\displaystyle \frac{D}{D_c}}\right)}^m$, Where ${\varphi }_{max}$is the upper-bound porosity at the critical fractal dimension$D_c$, and $m(~1-1.5)$ is an empirical sensitivity exponent calibrated from image-derived pore morphology.

This monotonic form captures the enhanced connectivity and reduced tortuosity associated with increasing $D$, consistent with the multiscale features observed in Section 2.3.2. The permeability-porosity relationship is formulated as a fractal scaling law [44,45].

Similarly,Permeability $k\left(x,t\right)$ follows a fractal-dependent scaling [44]:

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

(2.21)

Physically, this equation characterizes the hydraulic aspect of fractal behavior, where the exponent ${\displaystyle \frac{D\left(x,t\right)}{3-D\left(x,t\right)}}$captures the bidirectional feedback between fractal dynamics and flow capacity. Specifically, variations in $D\left(x,t\right)$modify both $\varphi \left(x,t\right)$ and $k\left(x,t\right)$,while stress-induced permeability alterations in turn reshape the evolving fractal geometry.

It is noted that $f_{\varphi }\left(D\right)$ in the empirical expression represents the geometric scaling of porosity with fractal dimension, whereas Eq. 2.21 reflects the hydraulic scaling between permeability and porosity. Together, these relations provide a physically consistent framework for linking pore-scale morphology and flow behavior within the DFDM model.

2.4.2. Boundary and Initial Conditions

For simplicity, the spatial evolution of $D\left(x,t\right)$ was modeled along the axial direction of the core plug (1D representation), which captures the dominant gradient of pore compaction and fluid transport.

Domain: The governing domain $\Omega =\left[0,L\right]$corresponds to a 1D spatial representation of the flow-through core plug,

Initial condition: The initial condition for the effective fractal dimension was prescribed as:${D\left(x,0\right)=D}_0+\xi \left(x\right)$, the effective fractal dimension field $D\left(x,t\right)$was initialized as a spatially varying baseline $D_0\left(x\right)=2.57$ obtained from the static grayscale maps, , and $\xi \left(x\right)$ is a small random perturbation $(<1\%)$representing natural micro-scale heterogeneity.

Boundary conditions: At the inlet $\left(x=0\right)$ and outlet $\left(x=L\right)$, Neumann-type boundary conditions: ${\displaystyle \frac{\partial D}{\partial x}}{|}_{x=0,L}=0$ was imposed, ensuring no net transfer of structural irregularity across the boundaries.

Stochastic term: The stochastic forcing term $\xi \left(x,t\right)$was implemented as a zero-mean Gaussian random field with short spatial and temporal correlation lengths $\left({\mathcal{l}}_c\ll L,{\mathcal{t}}_c\ll T\right)$, which emulates unresolved local perturbations due to microfracture growth and stress fluctuations. Here, ${\mathcal{l}}_c$and ${\mathcal{t}}_c$denote the spatial and temporal correlation scales of the stochastic forcing field, respectively, both much smaller than the characteristic domain length $\mathrm{L}$ and simulation duration $\mathrm{T}$, ensuring effectively uncorrelated (white-noise-like) perturbations. Representative values used in calibration satisfy ${\mathcal{l}}_c/L\approx 0.02~0.04$ and ${\mathcal{t}}_c\ll T\approx {10}^{-3}$, consistent with the image resolution and temporal discretization scales. The variance of the stochastic term was calibrated to match the experimental scatter of the measured fractal dimension ${\sigma }_D\approx 0.03$, as obtained from the static grayscale pore maps.

2.4.3. 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 ·\nabla Dd\Omega {=\int }_{\Omega }\omega \left(-\beta tr\left({\sigma }_{ij}\right)+\gamma ‖\nabla p‖+\xi \right)d\Omega$$

(2.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}$$

(2.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)}}$$

(2.23)

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 a quantitative link between experimental observation and DFDM model initialization. As illustrated in Figure 2, the workflow begins with thin-section micrographs of the rock matrix (a), which are segmented into binary pore-matrix representations (b). The corresponding gray-level surface is then reconstructed from the grayscale intensity field, where local variations encode pore-scale heterogeneity. This surface representation defines a pseudo-3D topological structure whose fractal dimension theoretically ranges between 2 and 3. The subsequent box-dimension regression using the MATLAB interface (c) yields a baseline fractal dimension of approximately $D_0\approx 2.57$, representing the surface roughness of the pore network. This baseline value is thus adopted to initialize DFDM predictions and serves as the experimental-theoretical bridge connecting static imaging with dynamic model evolution. 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}}$$

(2.25)

In general, empirical predictive models are expected to achieve$R^2>0.$9 to ensure statistical robustness. In contrast, the calibrated DFDM achieved $R^2=0.88$, which remains satisfactory for a physics-informed model constrained by first-principles dynamics and experimental uncertainties.

This result confirms that DFDM reliably reproduces both the structural constraints from thin-section analyses and the flow responses observed in 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. Figures Error! Reference source not found.-Error! Reference source not found. 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.

Model-Experiment Consistency of Fractal Dimension Thresholds 

The permeability-porosity relation expressed in Eq.2.21originates from the fractal scaling theory of porous media [44,45,46], where the exponent $D/\left(3-D\right)$characterizes the non-linear sensitivity of permeability to porosity as a function of the pore-network fractal dimension $D$.

Dependence of  $k$   on   $D$:

When $D$ increases toward 3, the exponent $D/\left(3-D\right)$ diverges, implying that small changes in porosity result in large variations in permeability. This reflects a highly complex, heterogeneous structure where flow becomes concentrated along preferential channels, corresponding to the mixed-phase formation stage in CO₂ flooding. Conversely, when D decreases toward 2, the exponent approaches unity, indicating a smoother and less connected pore system with unstable displacement and fingering (post-breakthrough stage).

Mathematical sensitivity analysis: 

The rate of change in the scaling exponent, $f^{\text{'}}\left(D\right)={\displaystyle \frac{3}{{\left(3-D\right)}^2}}$ ,increases sharply from $f^{\text{'}}\left(2.2\right)\approx 7$ to $f^{\text{'}}\left(2.6\right)\approx 19$, highlighting a distinct non-linear inflection zone between $D=2.2$ and $D=2.6$. Within this range, the coupling between porosity and permeability transitions from quasi-linear to strongly non-linear behavior. This window $\left(2.2<D\le 2.6\right)$ represents the mixed-phase displacement stage, characterized by effective connectivity and relatively uniform flow, i.e., the most stable sweep regime.

Quantitative and experimental consistency: 

The mathematical inflection region (2.2-2.6) coincides with the three displacement stages identified in the CO₂ flooding experiments (Figure 4). Specifically, $D>2.6$ corresponds to the mixed-phase formation stage with emerging preferential channels; $2.2<D\le 2.6$ corresponds to the mixed-phase displacement stage with uniform and stable sweep; and $D<2.2$ corresponds to the post-breakthrough stage with localized fingering.

Relation to  $D_0=2.57$:

The threshold $D_0=2.57$ marks the onset of non-uniform flow, where $f\left(D\right)=D/\left(3-D\right)\approx 6$. This indicates that permeability becomes approximately six times more sensitive to porosity variations than in the quasi-linear regime, signifying the initiation of channelization. Hence, $D_0$ lies within the mathematically derived and experimentally verified transitional window (2.2-2.6), providing quantitative justification for the adopted thresholds.

In summary, the D ranges correspond to physically distinct regimes of pore-network evolution and flow behavior, demonstrating the internal consistency between the DFDM formulation and the observed CO₂ flooding dynamics. Based on the mathematical sensitivity analysis of Eq.2.21 and the observed transition near $D_0=2.57$, the dynamic fractal dimension can be divided into three characteristic flow regimes, as summarized in Table 1.

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 (increasing D ) and flow-driven dissolution/enhancement (reducing 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.

Practical implications for subsurface engineering applications 

The distinction between compaction-dominated and flow-driven regimes also provides practical insight into subsurface engineering strategies. In CO₂ sequestration, compaction-dominated regimes generally enhance long-term storage security by reducing pore connectivity and permeability, thereby stabilizing trapped CO₂. Conversely, flow-driven regimes, while promoting dissolution and reactive transport, may increase leakage risks if not balanced by mechanical sealing. In enhanced oil recovery (EOR), moderate flow-driven behavior, manifested as controlled channelization, can improve sweep efficiency by creating preferential flow paths through low-permeability zones. However, excessive channelization or fingering leads to early breakthrough and reduced recovery. The DFDM framework thus provides a predictive tool to identify and manage the transition between desirable and detrimental flow regimes.

3.4. Limitations and Future Work

While DFDM demonstrates substantial predictive capability, several limitations merit consideration and suggest directions for further development:

• Computational cost:

The current DFDM simulations rely on high-resolution time-lapse imaging and fine spatial grids, resulting in considerable computational demand. This limits large-scale or long-term applications unless advanced solvers or parallelization strategies are employed.

• Representativeness of input data:

The model's predictive accuracy depends strongly on the representativeness of the input data. Reconstructions derived from 2D thin-section images may under-sample sub-resolution heterogeneities that are crucial for flow and transport, leading to potential scale-dependent biases.

• Parameter sensitivity:

The DFDM's sensitivity to stress, fracture-pore coupling coefficients, and boundary conditions can introduce uncertainties when extrapolating beyond laboratory settings. Further calibration under varied flow and stress conditions is necessary to improve robustness.

• Model extensibility:

The hierarchical and fractal-based architecture of DFDM allows natural integration with multiphase flow models, geochemical feedbacks, or physics-informed neural networks. Future efforts will focus on coupling DFDM with hybrid multi-scale solvers and machine-learning-assisted parameter estimation to enhance computational efficiency and generalizability.

Despite these limitations, DFDM provides a robust and physically grounded pathway to bridge pore-scale heterogeneity with reservoir-scale performance, offering a flexible foundation for future multi-physics and multi-scale simulations.

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 analysis, 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 DFDM predicts the proposed representative dynamic ranges of fractal dimension $D$ (2.2-2.6) were validated within the experimental ranges of effective stress (5-25 MPa), a constant injection rate of approximately 1.304 $ml·{min}^{-1}$, and system pressure varying from 30-42 MPa during CO₂ displacement. This configuration ensured steady-state flow and allowed pressure evolution to be monitored continuously along the core. The model reproduces the measured permeability-porosity-fractal trends with a mean error below 10% ($R^2=0.88$). These ranges define the model's current domain of applicability, extension to other fluids or mineralogies would require additional calibration.

The validated DFDM framework quantitatively captures stress-induced transitions in pore complexity during CO₂ flooding, linking stress-induced compaction and flow-driven reorganization with measurable changes in fractal dimension. The Model validation against laboratory experiments demonstrates strong agreement with dynamic 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 offers 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 scalable, quantitatively validated tool for energy recovery, carbon sequestration, and sustainable subsurface resource management. DFDM thereby unifies pore-scale imaging and reservoir-scale modeling, providing a robust framework to quantify and predict dynamic heterogeneity in subsurface energy systems.

Appendix A.1.Table of Symbols

Table A1. lists all symbols, parameters, and constants used in this study, together with their units and physical meanings.

Table A1. lists all symbols, parameters, and constants used in this study, together with their units and physical meanings.

Appendix A.2. Effective and Dynamic Fractal Dimension:Derivation and Implementation

Appendix A.2.1. Derivation of the Effective Fractal Dimension (Eq.2.6)

The pore radius distribution is assumed to follow a Weibull probability density function (PDF):

$$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]$$

where $m$is the Weibull modulus describing pore-size heterogeneity, and $r_0$ is the scale parameter characterizing the representative pore size.

To relate this statistical form to a fractal representation, we recall that the pore number density $N\left(r\right)$ in a fractal medium generally follows a power-law scaling:

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

Taking logarithmic differentiation of the probability density function gives a generalized definition of an effective fractal dimension as the local slope of $\mathit{log}p\left(r\right)$ with respect to $\mathit{log}r$:

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

Substituting the Weibull form into the above expression yields:

$$\mathit{log}p\left(r\right)=\mathit{log}m-\mathit{log}{r}_0+\left(m-1\right)log\left({\displaystyle \frac{r}{r_0}}\right)-{\left({\displaystyle \frac{r}{r_0}}\right)}^m$$

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

Thus,

$$D_{eff}\left(r\right)=m{\left({\displaystyle \frac{r}{r_0}}\right)}^m-\left(m-1\right)$$

For most natural pore systems, $m$ typically ranges between 1 and 5, and the exponential term dominates the variation with scale.

When the distribution is sufficiently broad $r\gtrsim r_0$,the offset $\left(m-1\right)$ is negligible compared to the exponential term, leading to the simplified expression adopted in the main text:

$$D_{eff}\left(r\right)=m{\left({\displaystyle \frac{r}{r_0}}\right)}^m$$

• Assumptions underlying this derivation

Self-similarity: the pore network exhibits scale-invariant statistics over the range where the Weibull form is valid.

Power-law number density: the pore population can be expressed as $N\left(r\right)\propto r^{-D}$, linking probabilistic and fractal descriptions.

Isotropy: pore orientations and shapes are statistically isotropic, allowing the use of a single scalar radius $r$.

Continuity: the distribution is smooth and differentiable with respect to $r$, enabling the use of local logarithmic derivatives.

• Physical interpretation

$D_{eff\left(r\right)}$ quantifies the local scaling exponent of the pore-size distribution on a log–log plot. In the asymptotic limit ${r\ll r}_0$, the Weibull distribution approaches $p\left(r\right)~r^{m-1}$, giving $D_{eff}\approx m$. As $r_0$ and $m$ evolve under stress and fluid–rock interaction, $D_{eff}$ becomes a dynamic descriptor reflecting the instantaneous geometric complexity of the pore network.

Appendix A.2.2. Parameter Extraction Algorithm and Stability Criteria

The parameter extraction procedure was developed to implement the Dynamic Fractal Dimension Model (DFDM) introduced in Section 2.3.1 and defined by Eqs. (2.12-2.14).

These equations describe the temporal evolution of the fractal dimension $D\left(t,x\right)$through a generalized scaling relation that links the measured pore-space occupation $N\left(\epsilon ,t\right)$to time-dependent structural parameters $C_0\left(t\right)$,$a$, and the fluctuation term $\eta \left(\epsilon ,t\right)$. The section provides the detailed algorithm for estimating the parameters, defining the stability criteria, and ensuring reproducible evaluation of $D\left(t\right)$ across temporal sequences.

The parameters $C_0\left(t\right)$,$a$, and $\eta \left(\epsilon ,t\right)$ are estimated from time-series pore images or segmentation maps via a multi-scale regression procedure, summarized below.

• Algorithm (pseudocode)

Input: Binarized pore images $I(x,t)$ across $N$ time frames

Output: $D\left(t\right),C_0\left(t\right),\alpha ,and\eta (\epsilon ,t)$

1. For each time frame $t$:

 a. Compute $N\left(\epsilon ,t\right)=$ number of occupied boxes for a set of scales $\epsilon \_i$ ($\epsilon \_i$ logarithmically spaced).

 b. Perform linear regression of log [$N\left(\epsilon ,t\right)$] vs.$log(1/\epsilon )$within a stable scale window [$\epsilon \_$min, $\epsilon \_$max].

 c. Obtain initial $D\left(t\right)=$ -slope, $C_0\left(t\right)=$ intercept.

2. Fit exponential decay $C_0\left(t\right)\approx C\_{init}^{*}exp(-\alpha t)$to estimate$\alpha$.

3. Compute residual $\eta (\epsilon ,t)={N\left(\epsilon ,t\right)-C}_0\left(t\right)exp(-\alpha t){\epsilon }^{\wedge }\left\{-D\left(t\right)\right\}$C0(t) exp(-α t).

4. Evaluate noise level $\sigma \_\eta =std\left(\eta \right)/mean\left(N\right)$; if $\sigma \_\eta >0.05$, readjust $\epsilon$-range by excluding small scales.

5. Smooth $D\left(t\right)$by Savitzky-Golay filter (window $=$3-5 frames) to suppress temporal noise.

6. Repeat across all t to construct $D(t,x)$field.

• Stability and Sensitivity Criteria

• Scale window selection: A minimum of one decade of $\epsilon$(i.e.,$log(\epsilon \_max/\epsilon \_min)\ge 1$) is required for robust regression. Windows are chosen by maximizing the coefficient of determination $R^2>0.98$.
• Temporal stability:The time step $\Delta t$is limited by the imaging frequency such that$\left|D(t+∆t)-D\left(t\right)\right|<0.05$to ensure quasi-stationary evolution.
• $\epsilon$-range sensitivity: The variation of $D$with respect to $\epsilon$is evaluated as

$S_{\epsilon }=\left|\partial D/\partial log\epsilon \right|<0.03$,

ensuring that the estimated $D$is scale-independent within the selected range.

These criteria provide reproducible and noise-resilient estimates of the dynamic fractal dimension and associated parameters, fully consistent with the DFDM theoretical framework presented in Section 2.

Appendix A.3. Wavelet and Box-Counting Implementation Details

The Morlet wavelet was employed as the mother function for the continuous wavelet transform, owing to its superior time-frequency localization and smooth oscillatory kernel suitable for multiscale pore-structure characterization. The analysis was conducted over a normalized scale range of $2^{-9}$ to $2^0$, corresponding to the smallest and largest box sizes (1/512-1/2 of the image dimension) used in the subsequent box-counting regression. This range ensures consistent spatial resolution across spectral and geometrical analyses. Mirror-symmetric boundary padding was applied to suppress edge artifacts, and local mean detrending was performed prior to spectral regression to ensure stationary residuals and numerical stability.

For the box-counting procedure, binary segmentation of the pore-matrix system was performed using adaptive grayscale thresholding (Otsu method, which automatically determined a normalized threshold of $T\approx 0.42$), followed by morphological filtering to remove isolated artifacts. A grid hierarchy of nine box sizes, ranging from 1/512 to 1/2 of the image dimension with power-law progression, was adopted. The logarithmic relationship between the box size and the number of occupied boxes was fitted using the least-squares method to obtain the baseline fractal dimension $D_0\approx 2.57$. These standardized preprocessing and regression protocols ensure the derived scaling parameters and their consistent initialization within the DFDM framework.

Appendix A.4. Numerical Discretization and Solver Configuration

Equation (2.19) was solved numerically using a finite-difference scheme on a uniform structured grid. The computational domain $\Omega$ was discretized into $∆x=∆y=1.0\times {10}^{-3}$m, with second-order central differences applied to spatial derivatives and an explicit forward-Euler time integration.

The time step $∆t=1.0\times {10}^{-4}$s was selected to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition：

$$∆t\le {\displaystyle \frac{{\left(∆x\right)}^2}{2{\kappa }_{max}}}$$

where ${\kappa }_{max}$denotes the upper bound of the fractal diffusivity. Under these parameters, the maximum Courant number remained $<0.4$, ensuring numerical stability. The system was iterated until the relative change in the global L₂ norm of $D(x,t)$ between consecutive time steps satisfied

$${\displaystyle \frac{{‖D^{n+1}-D^n‖}_2}{{‖D^n‖}_2}}<{10}^{-6}$$

Boundary conditions followed those introduced in Section 2.4.2. All simulations were implemented in MATLAB R2023a, with sparse-matrix operations and adaptive time-stepping checks to prevent overshoot or divergence. The stochastic term $\xi \left(x,t\right)$ was generated as a spatially correlated Gaussian noise field (zero mean,$\sigma =0.01$) using a spectral-synthesis method.

These configurations ensure numerical stability, convergence.