DIGITAL PETROGRAPHIC APPROACH TO CARBONATE RESERVOIR EVALUATION IN MATLAB

2.1. Study Area

The CCCP, operational since 1979 in Ukraine’s Lviv-Volyn coal basin, has served as a key facility for processing low-grade coal from surrounding mines. The substantial waste pile at CCCP covers approximately 0.8638 km² and rises to 68 m, encompassing about 48.89 Mm³ of waste at an average density of 2680 kg·m⁻³ (Bosak et al., 2020). The accumulated waste poses environmental risks because of its large volume and the possibility of pollutant leaching. The study area is in Western Polissya, a humid continental (Dfb) climate zone with cold winters and warm summers. Meteorological data from the Volodymyr weather station indicate significant annual precipitation variability of 511 - 829 mm (see Figure 1), which provides the range used later for precipitation-based projections.

Figure 1. Total annual rainfall observed at Volodymyr station (Ukraine) during 2018-2024 (URL 1).

The CCCP waste primarily consists of argillite, siltstone, sandstone, and coal-associated sediments. The predominant minerals-siderite (FeCO₃), pyrite (FeS₂), and chalcopyrite (FeCuS₂) - undergo oxidation, significantly contributing to environmental contamination through iron and sulphur mobilisation (Kochmar et al., 2024). As a consequence, these mineralogical processes promote the leaching of major ions and trace metals into surrounding hydrological systems (Karabyn and Kochmar, 2025). Leachate from the site may reach nearby surface waters such as the Western Bug River, which flows into the Vistula. As Vistula runs mostly through Poland, this transboundary connection requires coordinated monitoring and management.

Documented landslides, groundwater contamination, and episodes of spontaneous combustion in waste piles underscore the need for robust hazard-mitigation measures to protect regional ecosystems and communities (Kochmar et al., 2022; Petlovanyi et al., 2023). Understanding precipitation-driven leaching processes is essential for designing targeted environmental protection measures, improving waste management practices in Ukraine's coal basin.

2.2. Research methods

Column experiments and processing of results. To simulate the leaching dynamics of TDS from spoil rock, a laboratory-scale column experiment was conducted using an apparatus assembled at the Environmental Safety Laboratory, Lviv State University of Life Safety (Kochmar et al., 2024; Stepova et al., 2023). The test material was argillite collected from the CCCP spoil heap. Ten grab samples from different sectors of the heap were composited to improve spatial coverage, then air-dried at 20-25°C to constant mass over 30 d, gently comminuted in an agate mortar, and dry-sieved through a 1.25 mm stainless-steel mesh. The < 1.25 mm fraction was used in subsequent leaching tests, with no secondary fractionation. Sampling targeted 0.2-0.3 m to capture the biologically active and hydrologically responsive near-surface layer of the spoil, where infiltration, oxygen exchange, and temperature fluctuations are most pronounced, promoting sulphide oxidation, dissolution of readily soluble salts, and mobilisation of sulphate and chloride ions. A similar focus on near-surface layers has been used for coal-waste piles undergoing self-heating, where the upper horizons control infiltration and thereby increase the leaching of major ions and potentially toxic elements (Espinha Marques et al., 2024, Anghelescu and Diaconu, 2024).

For the tests, two material conditions were used: unburned argillite; and thermally altered argillite obtained by laboratory firing of the composite material in a muffle furnace for 4 h at temperatures up to 600°C. In this study, “thermally altered (burnt) argillite” refers to argillite that has undergone oxidation and partial mineralogical transformation through laboratory firing at 600°C, producing a reddish-brown material. The term “unaltered (unburnt) argillite” denotes material that has not experienced thermal alteration and retains the original grey, compact lithology typical of the parent rock. Each experiment used a 100 g subsample packed into a vertical column 10 cm high and 36 mm in diameter (cross-sectional area 0.001018 m²; packed volume ≈101.8 cm³) (see Figure 2). To ensure reliability, each treatment was run in triplicate: three identical columns were packed in parallel with the same mass of material (burnt or unburnt) and operated under identical flow and temperature conditions. Leachate was collected from each column at corresponding time intervals and analysed separately for TDS; reported values are means across the three parallel runs, which were used for the precipitation-response modelling. The corresponding standard deviations did not exceed the instrument’s measurement accuracy.

The concentrations of major ions in the filtrate were determined using standard chemical analysis methods at the Environmental Safety Laboratory at Lviv State University of Life Safety (certificate of attestation No. RL 091/21, dated 11/30/2021). In particular, acid titration to the appropriate end-point pH following Standard Methods was used (HCO₃⁻); argentometric titration with silver nitrate and potassium chromate as indicator, according to the Mohr method, was applied (Cl⁻). EDTA complexometric titration with murexide and Eriochrome Black T was performed (Ca²⁺ and Mg²⁺), with magnesium calculated by difference. Gravimetric determination as barium sulphate after precipitation with Ba²⁺ and subsequent ignition and drying of the BaSO₄ precipitate was used (SO₄²⁻). Sodium and potassium (Na⁺, K⁺) were estimated by the ionic charge balance method, based on equivalents of the directly determined major ions (APHA, 2017).

The leaching system comprised a vertical filtration column and a circulation pump that continuously drove deionized water through the packed bed. Deionized water (8.5 L) was recirculated at 300 mL·min⁻¹ for 24 h, yielding a cumulative passed volume of ≈432 L (≈50 flushing cycles). Aliquots of the filtrate were collected at fixed intervals (every 4 h over 24 h), with additional early-time sampling during the first 2 h of circulation (Kochmar et al., 2024).

Figure 2. Laboratory setup showing the column apparatus used for TDS leaching experiments.

Filtrate samples were collected at multiple time intervals and analysed for TDS (mg·L⁻¹), which were derived from measurements of electrical conductivity (EC). EC was measured with a bench-top EC 215 (Hanna Instruments), which automatically converts EC to TDS using a built-in conversion factor. The meter was two-point calibrated daily with standards bracketing the sample range and verified against a mid-range check, with recalibration performed if drift exceeded ±2%. The manufacturer-stated accuracy is ±1% of full scale. TDS values were calculated from temperature-corrected EC at 25°C according to the relationship TDS = f × EC₂₅, with f = 0.50 appropriate for sulphate-rich matrices. The method detection limit was 1.5 mg·L⁻¹ (based on 3σ of low-level EC checks at 25°C, converted to TDS). Quality control included reagent blanks, duplicate aliquots (~10%), and periodic check-standards (every 2 h), all within predefined acceptance limits.

Concentrations were multiplied by the cumulative filtrate volume to calculate cumulative mass (M) at each sampling stage; dividing M by the passed volume V yielded the mean leachate concentration (mg·L⁻¹). The results were used to model both the cumulative leaching behaviour and the incremental leachate concentration per unit water volume. The latter was obtained by calculating the change in TDS-derived mass over the change in water volume between successive sampling points, providing a dynamic view of concentration evolution throughout the experiment. Based on the observed M-V relationship, a linear regression model was fitted to describe salt removal as a function of water input. The model estimated M = a_v ·V + b_v , where M is cumulative leached mass and V is passed water volume. The regression coefficients a_v (mg·L⁻¹) and b_v (mg) were derived from experimental data and subsequently used to parameterise the site-scale precipitation model.

A two-stage approach was applied to process the column data. The primary experimental values were the TDS (mg·L⁻¹) in the leachate measured at specific time intervals. At the first stage, these values were converted into the cumulative mass of leached salts M (mg) by volume-weighting across successive sampling intervals:

[TeX:] \ M_{i}\ = \mathrm{\Delta}V_{i}{\sum_{j}^{}{TDS}}_{j}$M_i=\mathrm{\Delta }{V}_i{{\sum }_j^{}TDS}_j$ (1)

where ΔV_i is the incremental water volume between samplings.

At the second stage, the incremental leachate concentration between successive points was computed as:

[TeX:] \left( \frac{dM}{dV} \right)_{i} \approx \frac{M_{i} - M_{i - 1}}{V_{i} - V_{i - 1}}${\left(\frac{dM}{dV}\right)}_i\approx \frac{M_i-M_{i-1}}{V_i-V_{i-1}}$, (mg⋅L⁻¹) (2)

Calibration of column-derived regression coefficients for model parametrisation. Column data were used solely to calibrate the linear mass-volume relation M=a_vV+b_v for each material (burnt vs unburnt), yielding a_v (mg·L⁻¹) and b_v (mg). These dimensional coefficients are then transferred to the site scale in the mathematical model, where the real spoil area S_site, spoil mass m_ore,site, hydrologic partitioning (V_p,site,V_f,site) and annual duration T=1 yr govern the predicted mass export as a function of precipitation. To avoid ambiguity with the dimensionless formulation used later, we denote the column-regression coefficients by a_v, b_v here; when introduced into the site-scale dimensionless expression, they appear through dimensionless combinations denoted as a, b.

Column data were used solely to calibrate the linear mass-volume relation M=aV+b for each material (burnt vs unburnt), yielding the regression coefficients a_v (mg·L⁻¹) and b_v (mg). These coefficients are then transferred to the site scale in the mathematical model, where the real spoil area S_site, spoil mass m_{ore, site} hydrologic partitioning (V_{p, site,} V_{f, site}) and annual duration T=1 yr govern the predicted mass export as a function of precipitation. Thus, while the column analysis is expressed in terms of cumulative water throughput (M, V), the field application uses the calibrated (a, b) with site-specific geometry and hydroclimate.

Statistical processing of column data. For each treatment (unburnt and burnt), the cumulative mass–volume relation was fitted by ordinary least squares (OLS) as [TeX:] M = a_{v}V + b_{v}$M=a_{v}V+b_v$. We report the slope [TeX:] a_{v}$a_v$, intercept [TeX:] b_{v}$b_v$, their standard errors ([TeX:] SE_{a},SE_{b}$S{E}_a,S{E}_b$), the coefficient of determination ([TeX:] R^{2}$R^2$) and the sample size ([TeX:] n$n$). Two-sided 95% confidence intervals (CI) for the slope were computed as [TeX:] a_{v} \pm t_{0.975,\text{ }n - 2}\text{ }SE_{a}$a_v\pm t_{0.975,n-2}S{E}_a$, where [TeX:] t_{0.975,\text{ }n - 2}$t_{0.975,n-2}$is the Student’s [TeX:] t$t$-critical value for [TeX:] n - 2\ $n-2$degrees of freedom. Units are expressed in SI; since [TeX:] 1\text{ }mg\text{ }L^{- 1} \equiv 1\text{ }t\text{ }Mm^{- 3}$1mg{L}^{-1}\equiv 1tM{m}^{-3}$, the coefficient [TeX:] a_{v}\ $a_v$in mg·L⁻¹ is numerically identical to a per-volume sensitivity in t·Mm⁻³. Residual plots were inspected to check for obvious non-linearity or heteroscedasticity. The calibrated coefficients and their uncertainties are then propagated to per-100 mm sensitivities and site-wide annual increments using scalar error propagation (delta method), i.e. [TeX:] SE(K_{100}) = V_{100}\text{ }SE_{a}$SE\left(K_{100}\right)=V_{100}S{E}_a$and [TeX:] CI_{95\%}(K_{100}) = K_{100} \pm t_{0.975,\text{ }n - 2}\text{ }SE(K_{100})$C{I}_{95\%}\left(K_{100}\right)=K_{100}\pm t_{0.975,n-2}SE\left(K_{100}\right)$, where [TeX:] V_{100}$V_{100}$is the effective infiltration volume associated with +100 mm precipitation. These procedures follow standard treatments of simple linear regression and interval estimation (Montgomery et. al., 2021).

Mathematical Model Formulation. Building on column-leaching formulations and hydroclimatic transport-reaction metrics from prior work (Long et al., 2020; Kumar et al., 2020), together with similarity π-theorem scaling and leaching time-scale criteria (Vladyko et al., 2022; Crundwell, 2005), we formulate a dimensional-analysis framework for cumulative TDS leaching from coal spoil. We first outline the controlling processes and variables retained in the analysis, before assembling them into a compact set of dimensionless groups that couple precipitation forcing and infiltration with pore-scale retention/adsorption and heap geometry, thereby enabling scalable links from column conditions to field-scale predictions (Long et al., 2020; Kumar et al., 2020; Crundwell, 2005; Vladyko et al., 2022).

Leaching of SO₄²⁻, Cl⁻ and HCO₃⁻ from coal waste is controlled by meteoric infiltration, mineral structure, and oxidation at the surface. To generalise these effects, we applied dimensional analysis and similarity theory, producing a predictive model of cumulative ion release.

Under the working assumption that the cumulative mass of leached TDS scales with the infiltrated water volume and is modulated by physical characteristics of the spoil heap, we construct a system of dimensionless groups governing the leaching process. These comprise: (i) terms characterising the dependence of TDS mass on precipitation intensity and cumulative volume; (ii) groups combining hydraulic conductivity, leaching duration, porosity, density and water-solid interaction time; and (iii) geometrical and structural descriptors of the spoil that influence transport pathways and retention capacity.

The general non-dimensional form of the model is given by:

[TeX:] M_{TDS} = \ f\left( V_{p},\ V_{f},\ T,\ S,\ P_{p},\ m_{ore},\ \varepsilon,\ \rho,\ k_{ads} \right)$M_{TDS}=f(V_p,V_f,T,S,P_p,m_{ore},\epsilon ,\rho ,k_{ads})$ (3)

Where:

[TeX:] M_{TDS}$M_{TDS}$ – total mass of leached TDS over the integration period T within the chosen control volume (column or spoil block) (kg),

[TeX:] V_{p}$V_p$ – cumulative precipitation volume that falls onto the active infiltration area (S) during (T) (m³),

V_f – effective infiltration volume that actually reaches the reactive pore space (after run-off/evaporation losses) (m³),

T – leaching duration (integration time) (h or d),

S – infiltration area, i.e. the active surface of the spoil (or cross-section in a column set-up) through which water enters the system (m²),

[TeX:] P_{p}$P_p$ – precipitation intensity, (mm·h⁻¹),

[TeX:] m_{ore}$m_{ore}$ – mass of solid spoil within the control volume that participates in leaching over (T) (kg),

ε – effective (mobile) porosity, the fraction of pore space available to flow, dimensionless (0-1),

ρ – bulk density of the spoil mass, (kg·m⁻³),

[TeX:] k_{ads}\ $k_{ads}$– effective adsorption/distribution coefficient (m³·kg⁻¹).

The cumulative precipitation volume is defined as:

[TeX:] V_{p} = S\int_{0}^{T}{P(t)dt}$V_p=S{\int }_0^{T}P\left(t\right)dt$ (4)

where P(t) is the precipitation intensity (m·h⁻¹), S is the active infiltration area (m²), and T is the integration period (h or d).

For periods with approximately constant mean precipitation intensity Pₚ, this integral can be approximated as:

[TeX:] V_{p}\ \approx \ T \cdot S \cdot P_{p}\ $V_p\approx T\cdot S\cdot P_p$ (5)

The effective infiltration volume is then defined as V_f = η·Vₚ, where 0 ≤ η ≤ 1 is the infiltration coefficient accounting for runoff and evaporation losses:

[TeX:] V_{f}\ = \ \eta \cdot V_{p}\ $V_f=\eta \cdot V_p$ (6)

Assuming separable contributions from water flux and solid-phase reactivity, the normalised expression is:

[TeX:] \frac{M_{TDS}}{m_{ore}} = f(\frac{T\ \cdot \ S\ \cdot \ P_{p} \cdot \ m_{ore}}{V_{p}\ \cdot V_{f}\ \cdot \ \rho\ \ \cdot \ \varepsilon}) + \ \varphi\left( \frac{T\ \ \cdot \ S\ \ \cdot \ P_{p}\ \cdot \ m_{f}}{V_{p}\ \cdot V_{f}\ \cdot \ \rho\ \ \cdot \ \varepsilon} \right)$\frac{M_{TDS}}{m_{ore}}=f\left(\frac{T\cdot S\cdot P_p\cdot m_{ore}}{V_p\cdot V_f\cdot \rho \cdot \epsilon }\right)+\phi \left(\frac{T\cdot S\cdot P_p\cdot m_f}{V_p\cdot V_f\cdot \rho \cdot \epsilon }\right)$ (7)

Where:

[TeX:] m_{f}$m_f$ - mass of the sorbing (fine) fraction (kg).

The non-linear function [TeX:] \varphi( \cdot )$\phi \left(\cdot \right)$, associated with sorption and secondary mechanisms, can be expanded as:

[TeX:] \varphi( \cdot ) = \ c\ \cdot \frac{T\ \cdot \ \ S\ \ \cdot \ P_{p}\ \cdot \ m_{f} \cdot \ k_{ads}}{V_{p} \cdot \ V_{f} \cdot \ \rho\ \cdot \varepsilon\ } + \ d \cdot \frac{T\ \cdot \ S\ \cdot P_{p}}{V_{p} \cdot \ V_{f}} + \ e\ \cdot \ k_{ads}$\phi \left(\cdot \right)=c\cdot \frac{T\cdot S\cdot P_p\cdot m_f\cdot k_{ads}}{V_p\cdot V_f\cdot \rho \cdot \epsilon }+d\cdot \frac{T\cdot S\cdot P_p}{V_p\cdot V_f}+e\cdot k_{ads}$ (8)

Under the assumption that [TeX:] f(\frac{T\ \cdot \ S\ \cdot \ P_{p} \cdot \ m_{ore}}{V_{p}\ \cdot V_{f}\ \cdot \ \rho\ \ \cdot \ \varepsilon})$f\left(\frac{T\cdot S\cdot P_p\cdot m_{ore}}{V_p\cdot V_f\cdot \rho \cdot \epsilon }\right)$ is approximately linear in its argument, the linear component of the model (the first term of Equation 6) is:

[TeX:] \ a\ \cdot \frac{T\ \bullet \ S\ \bullet \ P_{p} \bullet \ m_{ore}}{V_{p} \bullet \ V_{f} \bullet \rho\ \bullet \varepsilon}\ \ \ \ $a\cdot \frac{T•S•P_p•m_{ore}}{V_p•V_f•\rho •\epsilon }$ (9)

where a, b, c, d, e are empirical coefficients obtained from laboratory column leaching experiments.

The function φ(·) is infinitesimal compared with the function f and may therefore be neglected. Therefore, the model can be simplified to the linear approximation:

[TeX:] \frac{M_{TDS}}{m_{ore}} \approx \ a\ \bullet \frac{T\ \bullet \ S\ \bullet \ P_{p} \bullet \ m_{ore}}{V_{p} \bullet \ V_{f} \bullet \ \rho\ \bullet \ \varepsilon\ } + \ b\ \ $\frac{M_{TDS}}{m_{ore}}\approx a•\frac{T•S•P_p•m_{ore}}{V_p•V_f•\rho •\epsilon }+b$ (10)

Whereas this formulation was validated primarily for TDS in our column experiments, the same framework is directly transferable to other major ions, provided their leaching behaviour is similarly governed by water input and physical retention. The resulting model enables estimation of contaminant fluxes (expressed, for instance, in kg·Mm⁻³, equivalent to 10⁻⁹ kg·m⁻³) based on precipitation regime, heap structure, and material characteristics, supporting applications in mine site risk assessment and remediation planning. These column‑scale metrics provide the empirical basis for the precipitation‑based scaling M_TDS = k·P , which is calibrated and applied to site-wide predictions.

3.1. General Leaching Patterns of Major Ions

This section reports laboratory outcomes and translates them into precipitation‑scaled predictions for site‑wide loads. A comparison of the concentrations of dissolved solids in the leachates from thermally altered (burnt) and unaltered samples shows a much higher concentration of TDS in the burnt material: after 24 h of leaching (8.5 L of deionised water), the TDS content reached 462 mg·L⁻¹ in the burnt sample, compared to 185 mg·L⁻¹ in the unburnt one (see Figure 3). In both materials, a pronounced early-time spike in TDS is observed at ~2 h, after which concentrations decline towards quasi-steady values over the remainder of the 24 h run.

This corresponds to more than double the cumulative mass of dissolved salts - approximately 3.93 g versus 1.57 g. These masses follow directly from the closed recirculating setup: M=∑C_jΔV_j over the 24-h run (see Equation 1), which is numerically identical to C_final×8.5 L because the entire dissolved load resides in the tank at the end; equivalently, M=(M/V)_avg×432 L (see Table 1).

Figure 3. Temporal leaching behaviour of TDS from argillite in a laboratory column test.

Differences in leaching behaviour are evident even at the early stages of the experiment. For instance, after only 2 h (≈0.7 L of water), the TDS in the burnt sample reached 183 mg·L⁻¹, which is 3.7 times higher than in the unburnt counterpart (50 mg·L⁻¹).

Analysis of individual ions in the final filtrates confirms this trend (see Figure 4). Sulphate (SO₄²⁻) concentrations in the burnt sample reached 318 mg·L⁻¹, more than double the value observed in the unburnt sample (141.5 mg·L⁻¹). Similar increases were observed for bicarbonates, Ca²⁺, Mg²⁺, K⁺, and Cl⁻.

Figure 4. Distribution of water-soluble ions in filtrates of thermally altered and unaltered argillite samples (sulphate bars scaled by ×0.1).

Conventions. Sulphate (SO₄²⁻) concentrations are shown rescaled by ×0.1 to maintain all ions on a common y-axis scale. This rescaling was applied solely for graphical clarity and does not affect the relative trends or ratios among ions.

3.2. Experimental Leaching as a Basis for Model Calibration

The laboratory column leaching test provided the empirical foundation for calibrating the mathematical model that describes TDS release from coal spoil material under rainfall conditions. Specifically, the experimental setup enabled computation of the average leaching of salt ions concentration per unit water volume (M/V, mg·L⁻¹) for unburnt and burnt argillite, as well as the cumulative leached mass. During the 24-hour experiment, 8.5 L of deionised water was passed through each 100 g sample under a controlled flow rate of 300  mL·min⁻¹. At the end of the flushing cycle, the cumulative TDS concentration in the filtrate reached 185  mg·L⁻¹ for the unburnt sample and 462  mg·L⁻¹ for the thermally altered (burnt) one. These values correspond to cumulative leached masses of approximately 1.57 g and 3.93 g per sample, respectively.

Assuming a linear relationship between water volume and salt release over the tested range, the average TDS per unit water volume was computed as M/V (units: mg·L⁻¹); detailed values are reported in Table 1. Given the geometric characteristics of the test columns (cross-sectional area ≈ 0.001 m²), the corresponding column-scale areal fluxes were approximately J_24h≈1570 g m⁻² day⁻¹ (unburnt) and J_24h≈3930 g m⁻² day⁻¹ (burnt), computed as J_24h=M/(A⋅24 h) using the column cross-section.

These experimentally derived column metrics (M, V, and M/V) informed the calibration of the precipitation‑based model, which was then used to simulate TDS leaching as a function of annual precipitation. Table 1 summarises the primary data and calculations used for model calibration. This ensured consistency between laboratory observations and the field-scale predictive framework.

Table 1. Column leaching results for unburnt and burnt argillite samples

(24 h circulation; 8.5 L deionised water; 100 g material).

Conventions. Average M/V is computed as total leached mass divided by the total passed volume (V_total = 432 L). Final TDS in filtrate refers to the concentration measured in the last collected solution portion.

To quantify linearity and propagate uncertainty into the site-scale model, the cumulative mass-volume relation was fitted as [TeX:] M = a_{v}V + b_{v}$M=a_{v}V+b_v$ for each treatment. The regression diagnostics-slope ([TeX:] a_{v}$a_v$), standard error of the slope ([TeX:] SE_{a}$S{E}_a$), 95% confidence interval, intercept ([TeX:] b_{v}$b_v$), its standard error, the coefficient of determination ([TeX:] R^{2}$R^2$) and sample size ([TeX:] n$n$) – are compiled in Table 2.

Table 2. Linear regression diagnostics for M–V calibration (column experiments)

Conventions. (i) Units: [TeX:] a_{v}$a_v$in mg·L⁻¹ are numerically equivalent to t·Mm⁻³ (since [TeX:] 1\text{ }mg\text{ }L^{- 1} = 1\text{ }t\text{ }Mm^{- 3}$1mg{L}^{-1}=1tM{m}^{-3}$). (ii) 95% CI computed as [TeX:] a_{v} \pm t_{0.975,\text{ }n - 2}\text{ }SE_{a}$a_v\pm t_{0.975,n-2}S{E}_a$(here [TeX:] t_{0.975,5} = 2.571$t_{0.975,5}=2.571$). (iii) Intercepts reflect linear extrapolation towards [TeX:] V \rightarrow 0$V\to 0$ and are not used for site-scale projections.

3.3. Precipitation-TDS response and sensitivity (per-volume and site-scale)

As a result of the mathematical modelling, we quantified the relationship between annual precipitation P (mm·yr⁻¹) and the annual TDS release per unit spoil volume M_TDS (t·Mm⁻³·yr⁻¹) for unburnt and burnt argillite (see Figure 5).

Figure 5. Sensitivity of per-volume TDS release to precipitation: comparison of unburnt vs. burnt argillite.

For clarity, the relations used in Figure 5 are written explicitly as:

M_{TDS unburnt}=3.31×10⁻³⋅P (11)

M_{TDS burnt}=8.78×10⁻³⋅P (12)

with slope parameters k_unburnt=3.31×10⁻³ and k_burnt=8.78×10⁻³. Across the tested precipitation range, the model yielded approximately 1.0-3.0 t·Mm⁻³·yr⁻¹ for unburnt and approximately 2.6-7.9 t·Mm⁻³·yr⁻¹ for burnt argillite, confirming a monotonic response and consistently higher release for thermally altered material.

In this precipitation-based representation (per-volume release vs. annual P), the burnt/unburnt slope ratio equals k_burnt/k_unburnt=2.653. At the column scale, however, the corresponding ratios obtained from independent metrics are close to ≈ 2.50: final filtrate concentrations were 462 vs. 185 mg·L⁻¹ (ratio 2.497), and cumulative leached masses were 3.93 vs. 1.57 g (ratio 2.503). Using a constant column cross-section of 0.001018 m², these masses imply areal fluxes of ≈ 3930 vs. ≈ 1570 g·m⁻² per 24 h (ratio 2.497). Taken together, the burnt/unburnt ratio is metric-dependent (concentration, mass, areal flux, or precipitation-based per-volume release) and should therefore be interpreted strictly within the normalisation used in each plot.

Model-based scaling using calibrated coefficients; areal annual fluxes Jyear are reported in the text for each scenario.

Site-wide scaling. Using the column-calibrated coefficients a, b and the site parameters (S_site, m_ore,site,,  ρ,  ε), the annual mass export of dissolved solids is estimated from precipitation as:

[TeX:] V_{p,site} = S_{site}P_{tot}$V_{p,site}=S_{site}{P}_{tot}$ (13)

[TeX:] V_{f,site} = {\eta V}_{p,site}$V_{f,site}={\eta V}_{p,site}$ (14)

where [TeX:] P_{tot}$P_{tot}$is the annual precipitation depth (m·yr⁻¹) and [TeX:] 0 < \eta < 1$0<\eta <1$ is the effective infiltration coefficient.

Results are expressed as mass load per year (e.g., t·yr⁻¹) and, where required, as areal annual flux:

[TeX:] J_{year} = \frac{M_{TDS,site}}{S_{site}\text{ }T}$J_{year}=\frac{M_{TDS,site}}{S_{site}T}$ (15)

(e.g. kg·m⁻²·yr⁻¹), enabling comparison across spoil heaps of different size and morphology. Formally, the areal annual flux is:

[TeX:] J_{year} = \frac{M_{TDS,\ site}}{S_{site}{\cdot 1}_{yr}}\ \ \ \ \ (kg \bullet m^{- 2}{yr}^{- 1})$J_{year}=\frac{M_{TDS,site}}{S_{site}{\cdot 1}_{yr}}(kg•m^{-2}{yr}^{-1})$ (16)

For S_site = 0.8638 km², the scenario values are: Dry (511 mm) – J_year ≈ 0.096 (unburnt) and 0.254 (burnt); Average (670 mm) – 0.1260 (unburnt) and 0.333 (burnt); Wet (829 mm) – 0.155 (unburnt) and 0.412 (burnt) kg·m⁻²·yr⁻¹.

For each additional 100 mm of annual precipitation, the predicted release increases by 0.331 t·Mm⁻³·yr⁻¹ for unburnt and 0.878 t·Mm⁻³·yr⁻¹ for burnt argillite. When scaled by the CCCP spoil volume V=48.89 Mm³, these per-100 mm sensitivities correspond to +16.18 t·yr⁻¹ (unburnt) and +42.93 t·yr⁻¹ (burnt) in site-wide annual loads. These sensitivities are constant within the applicability domain of the adopted linear approximation; slope values and their ratio will differ if a different metric or normalisation is used, which is why the specific definition is provided alongside each figure.

To operationalise the linear precipitation – TDS relations reported above, Table 3 summarises per-volume responses and the corresponding site-wide annual loads for three precipitation scenarios spanning the CCCP site range, while Table 4 lists the fixed marginal sensitivities per +100 mm of precipitation.

Table 3. Scenario TDS loads vs. annual precipitation (CCCP heap)

Conventions. (i) Units: per-volume values in t·Mm⁻³·yr⁻¹; site-wide totals assume V = 48.89 Mm³. (ii) Relations: M_TDS, unburnt = 3.31×10⁻³·P, M_TDS, burnt = 8.78×10⁻³·P; P in mm·yr⁻¹.

Table 3 summarises the per-volume model M_TDS = k·P and the corresponding site-wide annual loads for three representative precipitation scenarios – Dry (511 mm), Average (670 mm) and Wet (829 mm). These values confirm a monotonic increase with P and consistently higher magnitudes for burnt argillite. Relative to the Average scenario (670 mm), the Dry scenario reduces the site-wide load by − 25.73 t·yr⁻¹ (unburnt) and – 68.25 t·yr⁻¹ (burnt), whereas the Wet scenario increases it by +25.73 t·yr⁻¹ and +68.25 t·yr⁻¹, respectively. Across the full span from 511 to 829 mm, the site-wide TDS load rises by 51.46 t·yr⁻¹ (unburnt) and 136.50 t·yr⁻¹ (burnt).

Table 4. Fixed marginal sensitivity (linear model)

According to Table 4, the constant increments implied by the linear model per +100 mm of annual precipitation are summarised. This tabular presentation facilitates straightforward application of the model to varying climatic scenarios without the need for recalculation. These increments are invariant within the calibration domain and can be used to interpolate between scenarios or to adjust loads for intermediate P. They also facilitate quick updates using observed or forecast precipitation without re-fitting the model. For instance, +200-300 mm above mean annual precipitation would increase site‑wide TDS export by ≈32-64 t·yr⁻¹ for unburnt-dominant scenarios and ≈ 86-129 t·yr⁻¹ where burnt materials prevail.

Building directly on the per-volume precipitation - TDS relations (Equations 11-12; Figure 5), our results show a robust, approximately linear increase of annual salt release with increasing annual precipitation for both lithological states of argillite, with consistently higher magnitudes for thermally altered (burnt) material. The slope parameters k_unburnt = 3.31×10⁻³ and k_burnt = 8.78×10⁻³ (t·Mm⁻³·mm⁻¹) quantify this contrast and, together with the site‑scale volume of the CCCP spoil heap, allow direct translation from climatic variability to annual salt loads. This contrast between burnt and unburnt rock is consistent with the results of our earlier studies (Kochmar and Karabyn, 2023).

3.4. Ion-specific response of leaching to precipitation

Building on the linear precipitation-TDS relations described in Section 3.3, we apportioned the precipitation‑driven release to individual dissolved species using the final filtrate composition from the column experiments (see Table 5). For each material regime m (unburnt or burnt), we treat the leachate composition as quasi-stationary, i.e. c_i(t) ≈ α_i,m c_TDS(t). We compute α_i,m from the final tank composition as α_i,m=C_i,final/∑_j∈JC_j,final over the major ions listed in Table 5, so that ∑_iα_i,m=1. Integrating the mass balance M_i=∫c_i dV gives M_i=α_i,mM_TDS; together with M_TDS=k_TDS,mP and M_i=k_i,mP, this yields k_i,m=α_i,mk_TDS,m. For transparency, we also report ∑_iα_i,m and the variability of C_i/C_TDS across samples (median and coefficient of variation) to document the quasi-stationarity assumption. Here, α_i,m denotes the mass fraction of ion i among the major ions in the filtrate of material m (HCO₃⁻, SO₄²⁻, Cl⁻, Na⁺, K⁺, Mg²⁺, Ca²⁺), and k_TDS,m is the precipitation‑based slope (k_unburnt = 3.31 × 10⁻³; k_burnt = 8.78 × 10⁻³ t·Mm⁻³·mm⁻¹). The fractions are normalised to sum to unity, ensuring that Σ_ik_i,m = k_TDS,m. The resulting composition weights are: burnt argillite dominated by sulphate (≈79%) with bicarbonate ≈11%; unburnt argillite ≈72% sulphate and ≈13% bicarbonate; alkali and alkaline‑earth cations (Na⁺, K⁺, Mg²⁺, Ca²⁺) contribute only a few per cent each.

Table 5. Major-ion shares in the filtrate, used for partitioning the TDS slope

Conventions. (i) Fractions (α_i,m, %) are calculated from the major-ion concentrations and normalised within the major-ion set. (ii) For unburnt material, Cl⁻ is treated as 5.0 mg·L⁻¹ (upper bound).

To complement the composition weights in Table 5, Table 6 presents the marginal sensitivities per +100 mm of precipitation for individual ions, reported both as per-volume and site-wide values.

Table 6. Marginal increase in annual load per +100 mm of precipitation (ΔP = 100 mm).

Conventions.Values are reported both as per-volume (t·Mm⁻³·yr⁻¹ per +100 mm) and as site-wide increments (t·yr⁻¹ per +100 mm), assuming a spoil volume of V = 48.89 Mm³.

Building on these marginal sensitivities, Table 7 presents the predicted annual releases of major ions at a representative precipitation level of 670 mm·yr⁻¹, expressed both per-volume and as site-wide loads.

Table 7. Predicted major-ion releases at P = 670 mm·yr⁻¹

Conventions. Values are reported as per-volume rates ΔMᵢ/V (t·Mm⁻³·yr⁻¹) and as site-wide loads ΔMᵢ (t·yr⁻¹, in parentheses), computed for a spoil volume of V = 48.89 Mm³.

Sulphate dominates the precipitation-driven response in both materials, contributing ≈72% (unburnt) and ≈79% (burnt) of the bulk TDS slope. The bicarbonate fraction is 11-13%, consistent with carbonate buffering and dissolution. Among cations, Ca²⁺ and Mg²⁺ increments are modest but non-negligible. Since the model is linear, the difference in ion release between dry years (511 mm) and wet years (829 mm) can be obtained directly by multiplying the +100 mm sensitivities from Table 6 by the corresponding precipitation difference. As a reference point, Table 7 summarises the predicted annual releases of major ions at 670 mm·yr⁻¹ (the average precipitation scenario), showing that sulphate overwhelmingly dominates the total loads, followed by bicarbonate, with other cations contributing only marginally. This linearity assumption simplifies projections but may under-represent potential nonlinear responses during extreme hydrological events.

This equal-rate assumption follows the linear relationship M_TDS = kTDS, m·P. These methodological choices do not affect additivity or the dominant sulphate signal, but they should be reconsidered if time-resolved compositions become available. The results have direct implications for anticipating water‑quality impacts under extreme precipitation scenarios, which may escalate into local emergency situations requiring civil protection measures. Furthermore, the ion-resolved relationships established here underpin the scenario analysis linking climatic variability to composition-specific load changes. These quantified responses and sensitivities are interpreted in the context of mechanisms, published ranges, and implications for civil protection. These quantified responses and sensitivities are interpreted below in the context of mechanisms, published ranges, and implications for civil protection.

4.1. Leaching mechanisms

The approximately linear precipitation-TDS relation indicates that annual salt export is primarily water-flux-limited, whereas thermally altered argillite exhibits higher reactivity across scales. Two features are salient: (i) fixed marginal sensitivities enable scenario-based extrapolation across climatic regimes; and (ii) the magnitude of the burnt-unburnt contrast is contingent on the normalisation of the response metric (per spoil volume, per water volume, or per exposed area).

Thermal alteration augments salt release by breaking mineral lattices and exposing fresh reactive surfaces, by oxidising sulphides (e.g. pyrite) and promoting formation of readily soluble sulphate-bearing secondary phases (e.g. gypsum, anhydrite), and by altering pore architecture toward higher connectivity for percolating water (Castillo-Meza et al., 2020; Gao et al., 2025). Coal-spoil studies show higher SC/TDS ratios and sulphate-chloride fluxes from thermally altered (self-heated) materials in both column and field settings (Daniels et al., 2016b; Orndorff et al., 2015; Clark et al., 2017).

These pathways are consistent with the elevated sulphate and chloride measured in burnt filtrates and with widely observed increases in ionic strength from thermally affected mine wastes (Daniels et al., 2016a; Clark et al., 2017; Welch et al., 2021). In our columns, burnt material exhibits larger initial concentrations and cumulative masses than unburnt material under the same water input, supporting the interpretation that a thermally produced pool of labile salts is efficiently mobilised during wetting cycles. The persistence of a steeper precipitation-response slope in the annual model suggests that, beyond this transient pool, thermally altered rock maintains higher long-term reactivity. At other regional landfills, such as the Nadiya mine spoil tip, concentrations of major ions are also elevated, indicating that coal industry waste is a regional source of pollution throughout the Lviv-Volyn coal basin (Skrobala et al., 2022).

A sharp TDS peak around 2 hours (see Figure 3) reflects the rapid dissolution and wash-off of soluble secondary salts accumulated on reactive near-surface materials. The concentration then reduces as these salts are depleted and flow stabilises. Field and analogue studies of mine pools associate this ‘first flush’ to dissolution of efflorescent sulphate salts formed during oxidation. These salts produce highly mineralised waters at the onset of flushing, followed by a gradual reduction in concentrations (Mugova et al., 2024). Column tests on empty rock also show an ‘early weathering’ stage, characterised by high specific conductivity and peak ion concentrations, indicating rapid mobilisation of soluble salts prior to slower, reaction-limited release from less reactive phases (Martin & Langman, 2024). This interpretation aligns with the observed post-peak decline in our columns and supports treating early-time pulses as transient contributions superimposed on the longer-term, precipitation-modulated leaching captured by the linear model (Merritt & Power, 2022).

These findings align with laboratory and field evidence from other coal basins. Column studies on Appalachian spoils show high early-time TDS and specific conductance that decline toward quasi-steady levels after multiple pore-volume flushes, with higher values for unweathered, sulphide-rich mudstones and black shales (Orndorff et al., 2015; Daniels et al., 2016b). Work on closure salinity risk likewise documents rapid declines in Na-Cl where halite is initially present, followed by longer-term sulphate-dominated leaching (Park et al., 2013). Weathering-cycle simulations demonstrate that fine-grained fractions can produce extremely high initial TDS pulses that collapse after repeated wetting-drying cycles (Dang et al., 2019). At field scale, lysimeter and catchment monitoring confirm the same direction of change but over years to decades, with typical SC declining from >1000 - 3000 µS·cm⁻¹ to 300 - 700 µS·cm⁻¹, depending on lithology and management (EPA, 2011; Daniels et al., 2016a; Welch et al., 2021). Against this backdrop, our burnt/unburnt contrasts (≈2.5 at column scale; 2.653 at precipitation‑based per‑volume scale) fall within reported ranges.

4.2. Model adequacy and emergency risks

The linear precipitation-response model adopted in Section 3.3 aligns with established hydrological theory. In chemostatic catchments, solute concentrations vary only weakly with discharge, so annual loads scale approximately linearly with water flux; this behaviour is well documented and has been reaffirmed in recent reassessments (Godsey et al., 2009). At broader hydroclimatic scales, export vulnerability is governed by transport-reaction time-scale ratios (a Damköhler-type control), which implies stronger mass mobilisation with greater water inputs across realistic precipitation ranges (Kumar et al., 2020). These lines of evidence provide a physically consistent basis for a linear precipitation-mass relation and justify its application to the tested materials and site-relevant precipitation regime (Godsey et al., 2009, Knapp et al., 2024; Kumar et al., 2020).

From a civil protection perspective, the constant marginal sensitivities imply that wet‑year anomalies translate linearly into elevated TDS loads, with a stronger effect for burnt zones. Wet-year anomalies will proportionally increase site-wide loads, with stronger effects where burnt materials prevail. For instance, +200-300 mm above mean annual precipitation would increase site-wide TDS export by ≈32-64 t·yr⁻¹ for unburnt-dominant scenarios and ≈ 86-129 t·yr⁻¹ where burnt materials prevail. Since sulphate and chloride remain mobile under neutral to mildly alkaline conditions, such increases can influence ionic strength, cation exchange, metal speciation, and downstream aquatic ecology (Daniels et al., 2016a; EPA, 2011; Welch et al., 2021). These risks may escalate into local environmental emergencies, underscoring the need for preventive measures such as stabilisation of burnt sectors, interception of first-flush leachates, selective covers where f_b is high, and monitoring of SC and sulphate with thresholds informed by per-100 mm sensitivities.

4.3. Model limitations and future directions

The discussion reconciles column-scale evidence and a parsimonious precipitation-response model into a coherent picture of salt leaching from coal-derived argillites. Thermal alteration emerges as the principal intensifier, with effects that are visible across metrics and preserved after scale-up to annual, per-volume predictions. The linear form captures first-order climatic control while remaining transparent enough to support decision-making: every +100 mm of precipitation implies a fixed, material-specific increase in annual TDS export per-volume of spoil and, by extension, a predictable increment in site-wide load at the heap scale.

Methodologically, the approach integrates controlled column leaching (100 g samples, 24 h, 300 mL·min⁻¹ flow rate, ≈ 432 L cumulative throughput) with dimensional analysis and simple scaling to the spoil-heap volume (48.89 Mm³). The model thus links micro-scale leaching processes with meso- and macroscale hydrological forcing, allowing the evaluation of salinity export in a form directly applicable to early-warning and planning frameworks of civil protection.

The framework developed here is not limited to the Chervonohrad district. Its formulation is transferable to other spoil piles and mining waste facilities across Ukraine, particularly in eastern regions where spontaneous combustion, shelling, and extreme precipitation combine to create conditions of heightened risk. In such settings, thermal alteration and hydrologic forcing interact to intensify contaminant release, underscoring the importance of a predictive tool that can quantify expected load increments under different climatic scenarios. Extending application to eastern Ukraine provides an opportunity to support civil protection planning under conditions of ongoing military threat. The ability to forecast ion-specific releases strengthens preparedness for potential emergency situations, where contaminated runoff or groundwater could directly impact local communities and ecosystems.

Several limitations should be noted. First, depth coverage was restricted to 0.2 - 0.3 m, which targets the near-surface leaching zone rather than the full vertical variability of element mobility within the heap. This upper horizon is hydrologically responsive and biologically active - where infiltration, oxygen exchange and temperature fluctuations enhance sulphide oxidation, dissolution of readily soluble salts, and mobilisation of major ions (notably SO₄²⁻, Cl⁻, HCO₃⁻) into percolating waters - yet it does not resolve potential deeper pathways of mobilisation governed by compaction, moisture storage and redox stratification. Accordingly, the reported leachate chemistry should be interpreted as characteristic of the active near-surface source zone for first-flush export, with recognition that metre-scale subsurface processes - including low-temperature oxidation and self-heating fronts - may develop below the surface and alter migration pathways and fluxes at depth (Espinha Marques et al., 2024; Anghelescu & Diaconu, 2024).

The relationship between the sediment and the mass applied in this study represents a linear approximation of a more general model, in which non-linear sorption and secondary processes are grouped into a term [TeX:] \varphi( \cdot ).\ $\phi \left(\cdot \right).$In this analysis [TeX:] \varphi( \cdot )$\phi \left(\cdot \right)$ is set to zero, giving the simplified form used in further modelling. The approach was validated mainly for total dissolved solids using recirculating column experiments, and ion-specific loads were estimated by applying quasi-stationary composition weights. Additional uncertainty arises at field scale from hydrological representation and scale transition: runoff generation, preferential flow and seasonal storage can decouple precipitation from effective infiltration; stratification and patchiness within spoil heaps may bias scaling; and column geometries differ from field conditions. Accordingly, confidence intervals should be propagated when translating per-volume sensitivities into site-wide annual loads for risk management. Finally, the fixed per-100 mm sensitivities should be applied as constants only within the stated applicability and calibration domain of the linear model.

In summary, the synthesis highlights both the mechanistic insight and the applied relevance of the model. It provides a robust, adaptable framework for environmental management in mining regions and a practical foundation for integrated civil protection strategies across Ukraine and beyond.

Future work should aim to extend the modelling framework beyond the CCCP heap, applying it to other spoil piles and concentration plant dumps across Ukraine, particularly in the eastern regions that are subject to shelling and therefore present higher risks than the Chervonohrad district. Further developments should also integrate a dynamic water balance module that translates precipitation into effective infiltration as a function of cover, slope, and season. In order to implement such approaches under real constraints, optimisation based on dynamic programming is an effective method for selecting and sequencing technological actions (Lousada et al., 2024). In addition, model calibration should be extended not only to total TDS, but also separately to sulphates, bicarbonates and chlorides. Such improvements will enhance the accuracy and reliability of the model and will be useful in planning civil protection measures for mining regions.

4.4. Recommendations

Given the calibrated linear precipitation–TDS response (k_unburnt = 3.31 × 10⁻³ and k_burnt = 8.78 × 10⁻³ t·Mm⁻³·mm⁻¹) and the fixed marginal sensitivities per +100 mm of annual precipitation (+0.331 and +0.878 t·Mm⁻³·yr⁻¹; site-wide +16.18 and +42.93 t·yr⁻¹ for unburnt and burnt material, respectively), management should prioritise (i) stabilisation and hydrologic control of thermally altered (burnt) sectors, where per-volume slopes and loads are consistently higher; (ii) interception of first-flush leachates during heavy rainfall and snowmelt, as early-time pulses dominate short-term export; (iii) routine monitoring of specific conductance and sulphate as pragmatic early-warning indicators, with action thresholds derived directly from the +100 mm sensitivities in Table 4 and the precipitation scenarios in Table 3; (iv) selective covers and drainage measures in areas with high effective infiltration (η), to reduce Vf and thus annual loads; and (v) contingency planning for wet-year anomalies (+200 – 300 mm), which, under linear scaling, translate into predictable tens-of-tonnes increases in site-wide TDS export. These measures follow directly from the observed burnt/unburnt contrasts, the linear precipitation dependence over the site-relevant range, and the ion partitioning that highlights sulphate as the dominant contributor to annual loads.