Fractional Process with Jumps for Karstic Spring Discharge

1. Introduction

Karst spring discharge time series are known to exhibit highly irregular, heterogeneous, and multifractal dynamics [5,20]. These complex patterns reflect the underlying geologic heterogeneity, including the interaction between the porous matrix, fracture networks, and conduit systems [7]. Traditional linear hydrological models, which often rely on simplified assumptions of stationarity and Gaussian-like fluctuations, fail to capture the intricacies inherent in karst hydrology, such as long-range dependence, sudden jumps in flow rates, and the influence of extreme events [17].

The recession curves of karst spring hydrographs, frequently decomposed into distinct flow components (e.g., conduit-driven quick-flow and matrix-driven baseflow) serve as important indicators of aquifer properties, infiltration processes, and hydraulic connectivity [5]. Although conceptual models such as those proposed by Mangin [16] and Fiorillo [5] provide useful insights by distinguishing various storage and drainage compartments, they often lack the ability to represent the non-linear and fractal or multifractal nature of observed discharge records [16,20]. Empirical studies have shown that karstic systems display scale-invariant behavior: the complexity of flow paths and conduit geometries induces fractal patterns in spring discharge time series [23]. The fractal dimension D characterizes the degree of roughness and complexity on multiple scales [1,21], and thus serves as a crucial quantitative measure to understand karst processes [4,10].

Motivated by these findings, it is imperative to develop stochastic modeling frameworks capable of reproducing both the fractal scaling and the jump-like features observed in karst spring discharges. In this study, we employ a fractional Ornstein–Uhlenbeck process (fOUp) enriched with jump components to model and simulate karst spring discharge dynamics. The fOUp is a continuous-time Gaussian process with long-range dependence, characterized by a Hurst exponent H (with $0<H<1$), allowing it to capture persistent ($H>0.5$) and slowly decaying autocorrelations or anti-persistent ($H<0.5$) rough fluctuations inherent in karst systems [3].

The usual, mean-reverting Ornstein–Uhlenbeck process, (OUp) $X\left(t\right)$ is defined by the stochastic differential equation (SDE):

$$dX\left(t\right)=\kappa (\theta -X(t\left)\right)dt+\sigma dB\left(t\right),$$

(1)

where $\kappa >0$ is the mean-reversion parameter, $\theta$ is the long-term mean, $\sigma >0$ is the volatility, and $B\left(t\right)$ is a classic Brownian motion. Extending the OUp to a fractional context involves replacing $B\left(t\right)$ with a fractional Brownian motion $B_H\left(t\right)$ of Hurst index H. This introduces long-memory characteristics into the dynamics when $H>0.5$ and rough fluctuation when $H<0.5$.

The ultimate source for renewing the karstic hydrogeological system is infiltration from precipitation or from melting snow. To address sudden discharge spikes originating from infiltration, we incorporate a jump component into the driving force of the fractional Ornstein–Uhlenbeck (fOU) equation, resulting in a hybrid model that can be schematically written as:

$$dX\left(t\right)=\kappa (\theta -X\left(t\right))dt+\sigma d{B}_H\left(t\right)+dJ\left(t\right),$$

(2)

where $J\left(t\right)$ represents the jump component. In general modeling practice, often the jump component, $J\left(t\right)$, follows a compound Poisson process. However, it has the shortcoming that the jump times and jump sizes are independent, and the process has the Markov property, enforcing exponential interarrival times (i.e., times between jumps, also called as sojourn times in renewal theory). This property is counterintuitive in our application, since a long delay between jumps – induced mainly by long delay in precipitation events – leaves the aquifer drained. Hence, the water from the next precipitation event - or a quick snowmelt for its turn - is kept back in the aquifer, filling up the porous matrix and the fracture networks, allowing for only lesser peaks in spring water discharge. Inspired by renewal theory and supported by statistical evidence (see Section 6), we suggest using a semi-Markov jump process for $J\left(t\right)$ with Weibull-distributed interarrival times and also Weibull-distributed peaks [24]. Contrary to Markov models, semi-Markov models allow interdependence of the process’ state and the interarrival time, and essentially place no restriction on the interarrival time distribution, see, e.g. [24]. The proper correlation of the interarrival times and the jump sizes is secured by the probability integral transform and its inverse, along with the Gaussian copula. To the best of our knowledge, this model specification is novel in the literature. The jump component describes the abrupt short term, whereas the fractional component fine tunes the long-term behavior.

We apply this fractional jump Ornstein–Uhlenbeck framework to daily discharge data from karst springs in Northeast Hungary. The associated analysis illustrates the methodology’s capacity to replicate both the statistical and scaling features of the observed time series. By comparing simulated and observed discharge data, we demonstrate that this modeling strategy can reproduce essential empirical characteristics, such as fractal scaling, non-Gaussian tails, and sudden discharge extremes. Thus, the jump-fOUp approach provides new insights into the structure and complexity of karst aquifers and represents a promising avenue for understanding and predicting the non-linear, non-stationary, and fractal-like behavior of karstic spring discharges.

1.1. Outline of the Paper

The remainder of this article is organized as follows. Section 2 introduces the modelling framework—a fractional Ornstein–Uhlenbeck process augmented by a renewal–reward jump term—together with its mathematical properties and a short proof of well-posedness, which can be found in Appendix A. Section 3 reviews existing estimation techniques for the pure fOU process and adapts them into a two-tier procedure that later underpins our jump-fOUp inference. Section 4 describes the geological setting, the 19-year daily discharge dataset from four Hungarian karst springs, and the exploratory analyses that motivate our model choices. Section 5 presents the complete estimation–simulation workflow: detection of surges, copula-based generation of correlated jump pairs, fGn simulation for the fractional kernel, and the iterative scheme that couples both components. Section 6 compares observed and simulated hydrographs, fractal dimensions, and cumulative discharges, demonstrating that the jump-fOUp reproduces key statistical and visual features markedly better than short-memory alternatives. Section 7 discusses hydrological implications, the necessity of long-range dependence, and possible extensions such as time-varying parameters or precipitation-driven jump intensities. Finally, Section 8 concludes with the main findings and outlines directions for future research.

2. Fractional Ornstein–Uhlenbeck Process with Renewal Jumps

2.1. The Pure fOUp

The (mean reverting) fractional Ornstein–Uhlenbeck process (fOUp) extends the usual Ornstein–Uhlenbeck model by incorporating the fractional Brownian motion (fBm) as its driving noise term. To distinguish it from its jump-enriched counterpart, we refer to it as pure fOUp. The fractional Brownian motion generalizes the classic Brownian motion ($H=0.5$) to incorporate long-range dependence or rough paths. If $H>0.5$, the increments are positively correlated, leading to persistent behavior, i.e., long memory. If $H<0.5$, the increments are negatively correlated, resulting in rough paths and antipersistent behavior. Let $X\left(t\right)$ be a continuous-time stochastic process satisfying the stochastic differential equation (SDE)

$$dX\left(t\right)=\kappa (\theta -X\left(t\right))dt+\sigma d{B}_H\left(t\right),$$

(3)

where $\kappa >0$ is the mean-reversion rate, $\theta$ is the long-term mean toward which $X\left(t\right)$ reverts, $\sigma >0$ is the volatility parameter, $B_H\left(t\right)$ is a fractional Brownian motion with Hurst parameter $H\in (0,1)$. The solution of the fOU equation exists and is unique up to initial condition $X\left(0\right)$. It can be given by [3]:

$$X\left(t\right)={\mathrm{e}}^{-\kappa t}X\left(0\right)+\theta \left(1-{\mathrm{e}}^{-\kappa t}\right)+\sigma {\int }_0^t{\mathrm{e}}^{-\kappa (t-s)}d{B}_H\left(s\right).$$

(4)

For $H=0.5$, the fOUp reduces to the usual OUp driven by a classic Brownian motion. The long-memory or rough path behavior is inherited by the fOUp from fBm according to the Hurst parameter. Unlike Bm and fBm the fOUp has a stationary regime. In contrast to the usual OUp, the fOUp is non-Markovian for $H\ne 0.5$. For more details on the fOUp see, e.g., [3].

2.2. Adding Jumps to the fOUp

We extend the fOUp dynamics by including jumps in the driving force to capture the long-memory or rough behavior of hydrological processes alongside sporadic, strictly positive surges that can be extreme in magnitude. The fOUp component encodes persistent correlations or rough fluctuations observed in hydrology. The jump components are responsible for the shocks that arrive from extreme water infiltration. In the simplest analogue, the jump diffusion setting, one might add a compound Poisson process whose interarrival times are exponentially distributed and whose jump sizes can in principle have arbitrary distributions. However, in hydrological systems where only positive, heavy-tailed surges are relevant and where jump sizes and interarrival times can exhibit dependence, the assumption of a compound Poisson structure becomes too restrictive. Our approach uses an analogy to renewal theory and incorporates an additive renewal-reward term into the fOU equation.

We enrich fOUp with jumps by first specifying a renewal-reward process $J\left(t\right)$. The process $J\left(t\right)$ is specified by a sequence of increasing jump times $T_k={\sum }_{i=1}^k{\tau }_i,k=1,2,\dots$ with $T_0=0$. The interarrival times ${\tau }_i,i=1,2,\dots$ are positive random variables, and the jump sizes $S_i,i=1,2,\dots$ are positive random variables. We set $S_0=0$. With these $J(t,\omega )={\sum }_{i=0}^k{S}_i\left(\omega \right)$ if $T_k\left(\omega \right)<t\le T_{k+1}\left(\omega \right)$. Both interarrival times ${\tau }_i$ and jump sizes $S_i$ have a Weibull distribution in our application, and variables ${\tau }_i$ and $S_i$ are correlated. In this way, $J\left(t\right)$ is a point process, with a nondecreasing, right-continuous, step-function sample path structure, for which stochastic integrals can be defined in terms of standard Lebesgue-Stieltjes integration theory. As a result, its differential $\mathrm{d}J\left(t\right)$ can also be viewed as a stochastic measure that puts jump size weights $S_i$ on jump time points $T_i$. Now we define our jump-fOUp process $\tilde{X}\left(t\right)$ as:

$$\tilde{X}\left(t\right)={\mathrm{e}}^{-\kappa t}X\left(0\right)+\theta \left(1-{\mathrm{e}}^{-\kappa t}\right)+\sigma {\int }_0^t{\mathrm{e}}^{-\kappa (t-s)}d{B}_H\left(s\right)+{\int }_0^t{\mathrm{e}}^{-\kappa (t-s)}dJ\left(s\right).$$

(5)

Theorem 1. 

Let $\tilde{X}\left(t\right)$ be defined by equation (16). Then $\tilde{X}\left(t\right)$ satisfies the stochastic differential equation,

$$d\tilde{X}\left(t\right)=\kappa \left(\theta -\tilde{X}\left(t\right)\right)dt+\sigma d{B}_H\left(t\right)+dJ\left(t\right),$$

(6)

where, $\kappa >0$ is the mean reversion rate, θ is the equilibrium level, $\sigma >0$ is the volatility, $B_H\left(t\right)$ is a fractional Brownian motion (same as $W_H\left(t\right)$ used in appendix - standardized notation to $B_H\left(t\right)$), and $J\left(t\right)$ is the renewal-reward process defined above.

We leave the proof of the theorem in Appendix A.

The jump term $J\left(t\right)$ accumulates at random times governed by Weibull-distributed interarrivals and delivers jump magnitudes that are also Weibull-distributed. This construction ensures that each surge contributes a strictly positive, heavy-tailed (as the Weibull shape is invariably less than one in these applications, see Section 6) increment to $\tilde{X}\left(t\right)$, and that the negative correlation of interarrival times and jump sizes is preserved. The semi-Markov property emerges because the distribution of waiting times until the next jump depends on previously realized jump sizes (and vice versa), thereby departing from the memoryless exponential interarrivals of a compound Poisson process.

The flexibility of the model, resulting from the semi-Markov construction, allows the simultaneous reproduction of heavy-tailed surges and their further spreading effect and long-memory fluctuations. This resolves a common challenge in hydrological modeling, where only positive, subexponential jumps with significant extremal impact can occur. By merging Weibull-based jump processes with fractional diffusion, one captures rare yet sizeable excursions on top of a persistently correlated or anti-persistent baseline.

3. Parameter Estimation Methods

Much effort has been devoted to the parameter estimation problem for fOU processes. Assuming that a continuous record of observations is available in the non-mean-reverting case ($\theta =0$), Kleptsina and Le Breton [14] obtained first the formula for the drift parameter’s maximum likelihood (ML) estimation and proved almost sure convergence. Bercu, Courtin and Savy [2] proved a central limit theorem for the MLE in the case of $H>1/2$. They claimed without proof that the above convergence is also valid for $H<1/2$. Finally, Tanaka et al. [26] completed the case by studying the ML estimates when $\kappa >0$ and $H\in (0,1/2)$.

Turning to mean reversion with $\theta >0$ and keeping the assumption of continuous observation, Lohvinenko and Ralchenko [19] obtained the ML estimates of $\kappa$ and $\theta$ when $\kappa >0$ and $H\in (1/2,1)$.

Estimating the parameter vector $\Theta =(\kappa ,\theta ,\sigma ,H)$ of the fOUp in the more realistic case of discrete, time-equidistant observations ${\left\{X\left(t\right)\right\}}_{t=1}^N$ poses significant challenges due to the fractional nature of the driving noise. In the non-mean-reverting case specialized statistical methods were developed in [12,13], using the relationship between the Skorokhod (or divergence) and Stratonovich integrals. These works also provide limit theorems for the estimation as both the time horizon and the mesh size grow to infinity simultaneously, and the horizon increases faster than the mesh size. However, this condition is difficult to verify in a real application with finite sample size.

3.1. A Two–Tier Estimation Framework for the fOU Parameters

This section consolidates, in a single presentation, the continuous–record methodology of [12,13] and its exact, discrete–sample analogue due to [29]. Both procedures target the parameter vector $\Theta =(\kappa ,\theta ,\sigma ,H)$ of a fOUp

Tier I — continuous record ([13]).

Assuming full observation of the path ${\left\{X_t\right\}}_{0\le t\le T}$, the mean–reversion speed $\kappa$ can be consistently recovered by either least squares

$${\widehat{\kappa }}_{\mathrm{LSE}}=-{\displaystyle \frac{{\displaystyle {\int }_0^{T}X\left(t\right)dX\left(t\right)}}{{\displaystyle {\int }_0^{T}X{\left(t\right)}^{2}dt}}},$$

(7)

or the energy–type statistic;

$${\tilde{\kappa }}_{\mathrm{ET}}={\left({\displaystyle \frac{1}{{\sigma }^{2}H\Gamma \left(2H\right)T}}{\int }_0^{T}X{\left(t\right)}^{2}dt\right)}^{-1/\left(2H\right)},$$

(8)

where $\Gamma$ denotes Euler’s gamma function. Both estimators are strongly consistent and asymptotically normal as $T\to \infty$.

Tier II — exact discrete record ([29]).

Let $X\left(k\right)$ be an equally spaced sample of size N with the grid step $\delta =T/N,k=1,\dots ,N$. Define second-order differences and their quadratic variations.

$$D_{\delta }\left(k\right)=X(k+2)-2X(k+1)+X\left(k\right),S_2=\sum _{k=1}^{N-2}{D}_{\delta }{\left(k\right)}^2,S_1=\sum _{k=1}^{\lfloor (N-2)/2\rfloor }{D}_{2\delta }{\left(2k\right)}^2.$$

Closed-form, strongly consistent estimators are then obtained as

$$\begin{array}{cc}\hfill \widehat{H}& ={\displaystyle \frac{1}{2}}{log}_2\left(\frac{S_1}{S_2}\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \widehat{\sigma }& =\sqrt{{\displaystyle \frac{S_2}{(N-2)\left(4-2^{2\widehat{H}}\right){\delta }^{2\widehat{H}}}}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \widehat{\theta }& ={\displaystyle \frac{1}{N}}\sum _{k=1}^{N}X\left(k\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \widehat{\kappa }& ={\left({\displaystyle \frac{\frac{1}{N}{\sum }_{k=1}^N{\left(X\left(k\right)-\widehat{\theta }\right)}^2}{{\widehat{\sigma }}^2\widehat{H}\Gamma \left(2\widehat{H}\right)}}\right)}^{-1/\left(2\widehat{H}\right)}.\hfill \end{array}$$

(12)

4. Application to Karstic Spring Discharge Modeling

4.1. Geological and Hydrological Settings

The Aggtelek Karst, part of the Gömör-Torna Karst region on the Hungarian-Slovakian border, is predominantly composed of Triassic carbonates, mainly Wetterstein limestone and dolomites [27]. The Hungarian section comprises approximately 200 km² of east–west-oriented karst plateaus at elevations of 400–600 m a.s.l., dissected by valleys as low as 150–260 m. Karstification began in the Miocene under subtropical conditions, with significant Pliocene uplift that forms the present-day north-south topographic gradient [9,11]. The mean annual precipitation is about 600 mm, most of which infiltrates the plateau and resurfaces at karst springs. The area has a humid continental climate with a mean annual temperature of 9.1 °C, about 120–130 frost days, and 40 snow-covered days annually [27].

These geological and climatic data originate from Telbisz et al. [27]. The dataset, and descriptions therein are based on decades of multidisciplinary fieldwork and geological mapping [8,18], and are considered reliable and high-quality within the regional geoscientific community.

4.2. Data Description

The measurement of the daily spring discharge of the most significant karstic springs in the Aggtelek Karst was carried out by the Hungarian Water Resources Research Plc. (VITUKI) from the early 1960s to the 1990s [22]. We selected four springs (Csörgő, Kecskekút, Komlós, Lófej) with non-missing daily observations for 6200 days in the period from 1 March 1974 to 18 February 1993. The basic descriptive statistics of these decades of uninterrupted hydrological data series are shown in Table 1.

The springs are located within a relatively small area, close to each other (Figure 1). However, significant scale differences can be observed in their daily discharges: the multiannual averages vary between 287 and 14508.5 $m^3/\mathrm{day}$, see Table 1. Cave systems of varying lengths, are known to be connected with the springs [15,25]. In Jósva spring water comes to the surface from the 25 km long Baradla-Domica cave system. Komlós spring (Kom) is connected to Béke cave (approx. 7.2 km long). Szabadság cave (approx. 3.2 km) belongs to the aquifer of Kecskekút spring (Kec). Caves of much shorter confirmed lengths – typically a few hundred meters, as cave researchers put it to current knowledge - lie behind the Csörgő spring (Cso) [15]. The cave systems direct quick flows to the mouth of the springs, taking essential role in shaping the surges in the hydrographs.

5. Methodology: Parameter Estimation and Simulation

The objective is to model the daily discharge time series $X\left(t\right)$ of several karst springs. The model assumes that the discharge dynamics results from the superposition of a continuous, persistent background process and discrete, sudden jump events, often triggered by precipitation. A jump fOUp as described in Section 2 is employed. The analysis proceeds spring by spring; the following sections outline the methodology for a single representative time series $X\left(t\right)$, $t=1,\dots ,N$. The modeling process involves two main phases: estimating model parameters from the observed data and simulating the process using these parameters.

5.1. Phase 1: Parameter Estimation from Observed Data

5.1.1. Data Preparation and Initial Analysis

Let $X\left(t\right)$ denote the observed discharge at day t. The differenced series, $\Delta X\left(t\right)=X\left(t\right)-X(t-1)$ for $t=2,\dots ,N$, is calculated to highlight rapid changes or jumps (Algorithm A1, line 1).

5.1.2. Jump Component Identification and Characterization

• Jump Detection: Jumps are identified where $\Delta X\left(t\right)$ exceeds a high quantile threshold, $q_{\Delta X}=\mathrm{quantile}(\{\Delta X\left(t\right)\},p)$, where p is a high probability, chosen as 0.975 (Algorithm A1, line 2). Let $T_j$ be the time index of the j-th detected jump, $j=1,\dots ,M$ (Algorithm A1, line 3). (For heuristic explanation see Section 6.1).
• Jump Size Analysis: The jump sizes are $J_j=\Delta X\left(T_j\right)-q_{\Delta X}$ (Algorithm A1, line 4). Their distribution is analyzed by fitting candidate distributions (Weibull and Gamma) to $J_j$ using Maximum Likelihood Estimation (MLE) (Algorithm A1, lines 6). Goodness-of-fit is assessed using Kolmogorov-Smirnov (KS) tests. Let ${\widehat{k}}_J,{\widehat{\lambda }}_J$ denote the estimated jump size parameters.
• Jump Timing Analysis: The inter-arrival times ${\tau }_j=T_j-T_{j-1}$ (with $T_0=0$) are calculated (Algorithm A1, line 5). Their distribution is analyzed similarly, fitting a Weibull distribution to a slightly perturbed data, via MLE (Algorithm A1, line 8). Perturbation is necessary because being a continuous distribution, integers and ties do not fit the Weibull. Perturbation is created by adding a Uniform[0.01,0.5] variable to the jump sizes. The effect of the perturbation is then analyzed using 100 perturbations. Let ${\widehat{k}}_{\tau },{\widehat{\lambda }}_{\tau }$ denote the estimated jump time parameters.
• Jump Dependence Structure: The dependence between jump sizes $J_j$ and preceding inter-arrival times ${\tau }_j$ is quantified by the sample correlation coefficient $\widehat{\rho }=\mathrm{Corr}(J_j,{\tau }_j)$ (Algorithm A1, line 9). Typically, $\widehat{\rho }$ is significantly less than 0.

5.1.3. Background Process Characterization

• Hurst Parameter Estimation: The long-range dependence is quantified by the Hurst parameter $H\in (0,1)$. For the background process, assuming the fOUp structure, it is estimated from (9) (Algorithm A1, line 10). However, for the complete $X\left(t\right)$ process with jumps – when (9) is not applicable – it is estimated using the madogram method via the fractal dimension: $\widehat{H}=2-{\widehat{D}}_{madogram}$. It also serves to remain compatible with the findings of [4].
• Fractional Ornstein-Uhlenbeck (fOU) Parameter Estimation: The background process is modeled as an fOUp, and its parameters are estimated from a truncated series $X^{\prime }\left(t\right)$ (Algorithm A1, lines 11-12, see also Section 6).

  −

  Volatility ($\sigma$): Estimated from the second differences of $X^{\prime }\left(t\right)$ using $\widehat{H}$:

  $${\widehat{\sigma }}^2\approx {\displaystyle \frac{\mathrm{Var}\left({\Delta }^2{X}^{\prime }\left(t\right)\right)}{4-2^{2\widehat{H}}}}.$$

  (13)

  (Algorithm A1, line 13).

  −

  Mean Reversion Rate ($\kappa$): Estimated from the variance of $X^{\prime }\left(t\right)$ and $\widehat{\sigma }$:

  $$\widehat{\kappa }\approx {\left({\displaystyle \frac{\mathrm{Var}\left(X^{\prime }\left(t\right)\right)}{{\widehat{\sigma }}^2\widehat{H}\Gamma \left(2\widehat{H}\right)}}\right)}^{-1/\left(2\widehat{H}\right)}.$$

  (14)

  (Algorithm A1, line 14).

5.2. Phase 2: Simulation Setup

• Jump Simulation via Gaussian Copula: To preserve the correlation $\widehat{\rho }$ between jump times and jump sizes, a Gaussian copula is used. To create a large sample from the Gaussian copula, $L\gg M$ pairs of correlated standard normal variables $(Z_{1,i},Z_{2,i})$ with correlation ${\rho }_G$ are generated and transformed to uniform $(U_i,V_i)$ via the normal CDF $\Phi (·)$ (Algorithm A1, line 15-16). Simulated inter-arrival times ${\tau }_i^{\prime }$ and raw jump sizes $J_i^{\prime }$ are generated using the inverse CDFs of the fitted Weibull distributions (Algorithm A1, lines 20-29). Since this transformation does not preserve the correlation, ${\rho }_G$ is chosen slightly larger than $\widehat{\rho }$, the target one. Subsampling is then performed on these L pairs until we obtain a final sequence of $M_{sim}\approx M$ pairs $({\tau }_{sim,j},J_{sim,j})$, so many that jump times cover the [0,T] time interval, and their correlation reaches the target $\widehat{\rho }$ (Algorithm A1, lines 19-25). An innovation array $J^{*}\left(t\right)$ of length N is built by filling it up with the final jump increments $J_{sim,j}+q_{\Delta X}$ at the simulated jump times $T_{sim,j}={\sum }_{k=1}^j{\tau }_{sim,k}$, and with zeros elswhere. (Algorithm A1, lines 31-34).
• Fractional Gaussian Noise (fGn) Simulation: A sequence of fGn increments $ϵ\left(t\right)$ with Hurst parameter $\widehat{H}$ is generated by differencing a simulated fBm (e.g., by using the somebm package of R) and scaled using the estimated volatility $\widehat{\sigma }$ (Algorithm A1, lines 36-37).

5.3. Phase 3: Combined Process Simulation Loop

Writing $e^{-\kappa t}$ as $e^{-\kappa }·e^{-\kappa (t-1)}$ and integrating in (16) first from 0 to $t-1$ and then from $t-1$ to t we obtain

$$\tilde{X}\left(t\right)={\mathrm{e}}^{-\kappa }X(t-1)+\theta \left(1-{\mathrm{e}}^{-\kappa }\right)+\sigma {\int }_{t-1}^t{\mathrm{e}}^{-\kappa (t-s)}d{B}_H\left(s\right)+{\mathrm{e}}^{-\kappa t}{\int }_{t-1}^t{\mathrm{e}}^{\kappa s}dJ\left(s\right).$$

(15)

The last integral is $J\left(T_n\right)$ if $t=T_n$ for a jump time $T_n$, or zero otherwise. Since $\kappa$ is small in our case, we obtain an approximation

$$\tilde{X}\left(t\right)\approx {\mathrm{e}}^{-\kappa }X(t-1)+\theta \left(1-{\mathrm{e}}^{-\kappa }\right)+\sigma \Delta B_H\left(t\right)+J\left(T_n\right)·{\chi }_{\{t=Tn\}}.$$

(16)

Using the approximation the final discharge series $\tilde{X}\left(t\right)$ is generated iteratively for $t=1,\dots ,N$ :

$$\tilde{X}\left(t\right)=e^{-\widehat{\kappa }}·\tilde{X}(t-1)+\theta \left(1-{\mathrm{e}}^{-\kappa }\right)+\mathrm{Innov}\left(t\right)$$

(17)

where $\mathrm{Innov}\left(t\right)=\sigma ·ϵ\left(t\right)+J^{*}\left(t\right)$ is the total innovation ($\sigma ·$ fGn + jump), and $\tilde{X}\left(0\right)$ is the initial value (Algorithm A1, lines 38-42):

This procedure generates a synthetic time series $\tilde{X}\left(t\right)$ to replicate the statistical properties of the observed discharge.

6. Application: Jump Distribution Fitting and fOUp Parameter Estimation for Karstic Spring Discharges

6.1. Separating the Jumps

The surges in spring discharges are not merely isolated spikes; rather, they exert a dampening effect on the recession curve. Consequently, it is not sufficient to model these as simple jumps in the process itself — they must be incorporated into the driving force of the dynamics, specifically in the differential component of the stochastic differential equation (SDE). Therefore, to locate the jump positions and sizes, the differenced series has to be considered. The quantiles and maxima of the differenced spring discharge series are presented in Table 2.

The comparison of quantiles, maxima, and their relative proportions in Table 2 reveals the heavy-tailed nature of the underlying distribution. This observation rules out the adequacy of a pure fOUp for modeling. Replacing fractional Brownian motion (fBm) with a Lévy flight in the SDE could address heavy-tail behavior, but it would fail to capture the correlation between extreme values and their timing. This limitation motivated the development of our proposed approach.

To corroborate the jump threshold used, Table 2 lists the empirical $0.95$, $0.975$ and $0.99$ quantiles. Across all springs, the largest observed surge is roughly 20 to 450 times the $0.95$-quantile and 10 to 200 times the $0.975$-quantile, while the $0.975$-quantile is only 2 to 3 times the $0.95$-quantile. So, choosing the $0.975$-quantile level as the cut-off effectively separates the extreme outliers, interpreted as hydrogeological ’jumps’, from the long-memory background. 155 values are observed above the $0.975$-quantile level, and it is still suitable to fit a distribution to them with sufficient accuracy. Hence, we choose this threshold to determine the jumps.

6.2. Jump Distributions and Waiting–Time Dependence

We analyze karst spring discharges through two complementary layers. First, we model the jump component by semi-Markov dynamics, whose inter-arrival times and magnitudes follow Weibull laws. Second, after excising these jumps, the residual series is represented by a fractional Ornstein–Uhlenbeck baseline whose parameters are inferred from second–order differences as set out in Section 3. All statistical decisions (choice of distribution, selection of the jump threshold, and verification of long–memory parameters) are supported by formal goodness-of-fit tests and by the scaling diagnostics summarized below.

6.3. fOU Background: Continuous Versus Truncated Records

After removing every surge that exceeds the $97.5^{\mathrm{th}}$ empirical quantile, a fractional–Ornstein–Uhlenbeck (fOU) model is fitted to the resulting jump-filtered series. Table 3 lists the estimated volatility $\widehat{\sigma }$, damping coefficient $\widehat{\kappa }$, memory scale $e^{-\widehat{\kappa }}$ and Hurst exponent $\widehat{H}$.

Across all springs $\widehat{\kappa }$ is small, so the background discharge reverts slowly to its mean; the associated memory $e^{-\widehat{\kappa }}$$>0.93$ indicates pronounced persistence, while the conditional variance is far lower than in the raw data. The Hurst exponents fall between $0.675$ and $0.907$, confirming long-range dependence in the smooth baseline component.

6.4. Scaling Check via Fractal Dimension

To verify that the combined model—semi-Markov jumps plus fOU baseline—captures the multiscale roughness of the data, the fractal dimension D of each observed record is compared with the mean D from 100 simulated paths (Table 4). The close match confirms that the model reproduces the measured texture across scales.

The $H=2-D$ values in Table 4 describe the full jump + fOU process, so they are systematically lower than the background $\widehat{H}$ in Table 3: high-frequency jump spikes shorten the effective memory, whereas the jump-filtered baseline retains the longer-range dependence quantified in the first subsection.

Figure 2 displays the simulated and observed hydrographs of the four springs and demonstrates that simulations, based on the fitted model, exhibit strong visual agreement with the observed time series. This directly addresses a common criticism from hydrological and hydrogeological experts, who often contend that mathematical models—despite achieving good statistical fits in certain aspects—fail to produce realistic-looking hydrographs.

From a water resources management perspective, cumulative discharge is of paramount importance, as it directly reflects the volume of water available for various uses over a given period. Therefore, its accurate estimation is a key objective in hydrological modeling. Figure 3 demonstrates the performance of our model in this regard. To highlight the critical role of long-range dependence, we replace the background flow component with an autoregressive (AR) process, while retaining the same jump sequence used in the jump-fOUp simulation. The outcome is striking: the 100 cumulative discharge trajectories generated by the AR model (blue lines) significantly overshoot the observed discharge (red line), in some cases more than doubling it by the end of the period. In contrast, simulations based on the fOUp model closely track the empirical data, underscoring the model’s accuracy. Moreover, the variability among the jump-AR simulations is substantially higher than that of the jump-fOUp simulations, reflecting the latter’s smoother trajectories. These results illustrate that the jump-fOUp model yields more reliable cumulative discharge estimates, enhancing its utility for water resource planning and management.

7. Discussion

The jump-fOUp model integrates within a single stochastic framework the two primary controls of karst discharge: the long-range dependent background flow—comprising multiple hydrogeological components—captured via the fractional Ornstein-Uhlenbeck process, and the heavy-tailed recharge surges, modeled through a renewal-reward jump process with correlated, heavy-tailed inter-arrival times and magnitudes.

Figure 2 illustrates that simulated and observed hydrographs exhibit strong visual agreement, while Table 4 confirms that the model replicates the observed fractal dimension within the limits of sampling error. This demonstrates that the jump-fOUp model provides a generative mechanism capable of faithfully reproducing these key features.

In comparison to short-memory AR baselines [16], the jump-fOUp model offers significantly improved tracking of cumulative discharge volumes (Figure 3). This highlights that omitting long-range dependence—governed by the Hurst parameter H— can lead to potentially flawed assessments of water resources. When the background flow, represented by the truncated process, is instead modeled with a short-memory AR(1) process, a single AR parameter must account for all temporal dependencies. This results in systematically higher AR parameter estimates compared to the memory parameter in the fOUp model. Moreover, because this parameter acts not only on the driving noise but also on the jumps, it dampens the jump effects more slowly than in the fOUp case, thereby forecasting elevated discharges in the following days. The cumulative result—total discharged water volume—becomes unrealistically high, yielding an overly optimistic outlook for water resource management.

The improved accuracy in cumulative discharge estimation provided by our model arises from two key components:

(i)

a fractional kernel that propagates hydrological memory across months to years, and

(ii)

Weibull–Weibull jump pairs, whose negative correlation between size and waiting time reflects the natural depletion–recharge dynamics of the aquifer (Table 5).

Unlike Poisson jump-diffusions [28] or Lévy flights, the semi-Markov specification preserves both the heavy tail and the empirically observed dependence between surge size and lead time, a combination that classical multifractal surrogates describe but cannot generate mechanistically.

Parameter identification remains demanding, and further research into (e.g., iterative, likelihood-based, or Bayesian approaches for discretely observed jump-fOU processes would be beneficial. However, the current method yields parameters that produce simulations consistent with key empirical observations, including the crucial fractal dimension. Finally, extensions like time-varying $\kappa$ or precipitation-driven jump intensities—offer an avenue to embed climate signals directly into the discharge generator.

8. Conclusions

In the presented analysis a jump-fOUp with Weibull–Weibull semi-Markov surges reproduces the amplitude distribution, fractal dimension and cumulative discharge of four Hungarian karst springs. Long-range dependence ($H>0.5$) is indispensable: replacing the fractional kernel with short-memory noise inflates the accumulated volume by up to 50 %. Surge size and waiting time are negatively correlated; modelling them jointly is crucial for realistic flood peaks.

The framework therefore bridges descriptive multifractal analysis and process-based simulation, and can be extended with time-varying parameters or rainfall covariates for predictive use.