Comparative Studies of Nonlinear Models and Their Applications to Magmatic Evolution and Crustal Growth of the Huai’an Terrane in the North China Craton

Abstract

Power-law, inverse exponential and logarithmic models are widely used as empirical tools to describe anomalies in spatial and temporal geodynamic processes. However, the lack of clear interpretation of the relationships and distinctions among these models often makes their selection challenging, leaving them as empirical tools to be validated by data. This paper introduces these nonlinear functions derived from a unified differential equation, with parameters that reflect their relative nonlinearities and singularities, enabling their comparative application. By applying these functions to analyze magmatic events of the Huai’an Terrane, this study reveals two major crustal growth and reworking events between 2.6 and 1.7 Ga, each exhibiting distinctive nonlinear characteristics. The power-law function highlights strong nonlinearity and singularity during phases of intense magmatic activity, while logarithmic and exponential functions effectively characterize transitions between different tectonic processes. Geochemical data, including U-Pb zircon dating and Lu-Hf isotopic analyses, further validate the models by delineating distinct phases of crustal growth and reworking within the Trans-North China Orogen. The findings help connect the anomalies of frequency of magmatic events with the tectonic processes, providing important insights into the evolution processes of the North China Craton.

1. Introduction

In geoscience, anomalies can arise in various forms, such as unusual temperature patterns, flare-up magmatic activities, or seismic activities. Nonlinear processes of Earth systems, like tectonic plate interactions, volcanic activity and magmatic processes, contribute to these anomalies. Understanding these geology anomaly events is essential because they can reveal key insights into Earth’s dynamic processes, such as mantle convection, crustal deformation and the effects of plate subduction or collisions. Whether linked to tectonic movements, volcanic activity or subsurface composition, anomalies caused by nonlinear processes provide valuable insights that drive scientific discovery. The concept of local singularity, grounded in multifractal theory and fractal derivative theory, offers a powerful framework to describe the nonlinear properties of extreme geological events and associated anomalies [1]. Singularity in geological processes arises when small changes or inputs trigger disproportionately large effects, often due to nonlinear dynamics. Key mechanisms that cause such singularity in the Earth system include but are not limited to Multiplicative Cascade Processes (MCPs), Phase Transitions (PTs) and Self-Organized Criticality (SOC). MCPs such as diffusion, convection and turbulent flow in fluids or magma create feedback loops that amplify localized anomalies [2]. PTs involve change in state (e.g., solid–liquid–gas) that introduce sharp gradients and abrupt system changes. SOC includes systems like fault networks or sedimentary layers that evolve toward critical states, where minor triggers cause large-scale events. These anomalies are challenging to characterize through traditional statistical models due to their irregular, scale-dependent behavior and tendency for abrupt intensity changes across scales. Nonlinear models are essential for accurately capturing the inherent nonlinear patterns in these anomalies. This perspective enables a deeper understanding of the spatial and temporal complexities of extreme geological events, offering valuable insights into their behavior, prediction and management.

Indeed, one of the important properties of anomalies, especially in the context of geoscience, is the decay rate or pattern that describes how the anomaly diminishes or changes from the center toward the edge. This decay pattern often reveals key information about the underlying processes and mechanisms driving the anomaly. Understanding how an anomaly decays is crucial in interpreting its significance and potential impact. Several models or functions have been utilized to describe the decay regulations of anomalies. These include linear, logarithmic, exponential and power-law functions. However, the lack of clear interpretation of the relationships and distinctions among these models often makes their selection challenging, leaving them as empirical tools to be validated by data [3,4,5]. This gap can make it difficult for practitioners to understand how these functions are related and how they might be chosen for real-world scenarios and often based on validation of the datasets [6,7,8]. For example, it was stated in [8] that the power-law model remains one of the best-performing models, and other models that performed equally as well included the log-linear functional form in modeling concentration and discharge relations in hydrological modeling. The current paper demonstrates the derivation of these functions from a unified first-order differential equation, with parameters that highlight their distinct characteristics, including nonlinearity and singularity. Practical examples are provided to illustrate their applications to geological processes, with a focus on extreme magmatic events.

Different nonlinear functions such as inverse exponential, logarithmic and power-law models are applied to model the varying intensity and frequency of the magmatic events in the Huai’an terrane. Previous studies of tonalite–trondhjemite–granodiorite (TTG) gneisses and granites in the Huai’an terrane through U-Pb zircon dating, Lu-Hf isotopes and whole-rock geochemistry reveals a complex tectonic evolution, including periods of crustal growth (2.6–2.2 Ga), transition to crustal reworking (2.2–1.9 Ga) and crustal reformation (1.9–1.7 Ga) [9]. In the current paper, we demonstrate that in these phases, different nonlinear functions can be used to model the varying intensity and frequency of magmatic events. The results show that a power-law function can characterize the strong magmatic activity between 2.6 and 2.2 Ga, indicating a scaling relationship between event frequency density around a peak and time interval, with a distinct singularity. Logarithmic functions and exponential functions are used to analyze the transition phases where magmatic activity evolves from crustal growth to reworking, reflecting a gradual shift in the tectonic conditions. These models help connect the anomalies of frequency of magmatic events with the tectonic processes, providing important insights into the formation and evolution processes of the North China Craton (NCC).

2. Methods

In the context of nonlinear dynamics and fractal calculus, the density of an attribute within a dynamic system can be studied as a temporal process, with density variations over time. For an attribute with density $\mathsf{\rho }\left(t\right)$ at time t, this density is the result of underlying dynamic processes acting over time and can be modeled as $<\mathsf{\rho }\left(t\right)>$, which represents the expected density of the attribute at any given moment.

In a simple, linear temporal system where changes in density are solely due to uniform decay over time, we might assume a first-order decay model where the rate of change of $<\mathsf{\rho }\left(t\right)>$ with respect to time t is proportional to the density itself with a constant decay rate λ [10]:

$${\displaystyle \frac{d<\rho (t)>}{dt}}=-\lambda <\rho (t)>.$$

(1)

Solving the above equation yields the following exponential decay model:

$$<\rho (t)>=c{e}^{-\lambda t}.$$

(2)

This solution indicates that in a homogeneous, uniform system, attribute density decreases exponentially over time. Here, the system is smooth and linear, meaning that higher-order derivatives of $<\mathsf{\rho }\left(t\right)>$ will remain proportional to the original function, such that each derivative simply scales the density by (−λ)ⁿ, where n is the order of the derivatives. This temporal model is commonly observed in systems with uniform decay, such as radioactive decay in rocks or gravitationally driven water flow in a simple drainage network [11].

However, when considering systems in heterogeneous media or influenced by complex cascade processes over time—such as intermittent convective pulses, variable diffusion rates, or cyclic fluxes—the decay behavior deviates from uniformity. In these cases, the temporal change rate λ itself may fluctuate, influenced by irregular or nonstationary processes that lead to faster or slower decay rates. For example, the change rate of density is proportional to the density itself $<\mathsf{\rho }\left(t\right)>$ with a variable decay rate (–λ(t)) that changes with a change in scale (t) [12].

$${\displaystyle \frac{d<\rho (t)>}{dt}}=-{\displaystyle \frac{1}{\lambda (t)}}<\rho (t){>}^{\gamma },$$

(3)

where λ(t) is the coefficient as a function of scale t and γ is an integer index of nonlinearity. To solve Equation (3), one needs to know the functional form of λ(t). It is usually not possible to know the exact expression of this function since it may have different forms depending on the complexity of the system. However, the coefficient function 1/λ(t) can be approximated using the Taylor series expansion of 1/[λ(1/t)] as follows [12]:

$${\displaystyle \frac{1}{\lambda (t)}}={\lambda }_0+{\lambda }_1{\displaystyle \frac{1}{t}}+{\lambda }_2{\displaystyle \frac{1}{t^2}}+o({\displaystyle \frac{1}{t^2}}),$$

(4)

where the components λ₀, λ₁ and λ₂ correspond to the coefficients of the first three orders of 1/t. The remainder term states that the error for the approximation approaches zero when t decreases to zero and the order of the error is smaller than the order of (1/t²). Using this approximation, one can express the differential model as follows:

$${\displaystyle \frac{d<\rho (t)>}{dt}}\approx -({\lambda }_0+{\lambda }_1{\displaystyle \frac{1}{t}}+{\lambda }_2{\displaystyle \frac{1}{t^2}})<\rho (t){>}^{\gamma }.$$

(5)

Further simplifying the model by setting γ = 1, the new model (5) represents the combination of three different types of forms: normal exponential, power-law and exponential of inverse time:

$$<\rho (t)>=c{t}^{-{\lambda }_1}{e}^{-{\lambda }_{0}t+{\lambda }_2{\displaystyle \frac{1}{t}}}.$$

(6)

Model (6) gives a new four-parameter model with the parameter c as constant and λ₀, λ₁ and λ₂ representing the coefficients of the first three increment terms, respectively. In model (6), both the first term, $e^{-{\lambda }_{0}t}$, and the third term, $e^{{\lambda }_2{\displaystyle \frac{1}{t}}}$, are monotonically decreasing exponential functions of t, whereas the second term $t^{-{\lambda }_1}$ represents the monotonically decreasing power-law function of t. This new model can be converted into a simple form by applying logarithmic transformation to both sides of (6):

$$Log\rho (t) = c-{\lambda }_{0}t-{\lambda }_{1}Log(t)+{\lambda }_2{t}^{-1}.$$

(7)

The above model indicates that a general system can be approximated by combining the three models with four parameters. The first exponential model characterizes the non-singular decay rate, whereas the other two models characterize the singular decay rate.

If we set γ = 0, the solution of model (5) results in a combination of three different forms: linear, logarithmic and inverse time:

$$<\rho (t)>=c-{\lambda }_{0}t-{\lambda }_{1}Log(t)-{\lambda }_2{t}^{-1}.$$

(8)

The primary difference between model (8) and model (7) lies in the logarithmic transformation of the density in model (7). These two solution forms, (7) and (8), indicate that the forms in solution (7) generally exhibit higher nonlinearity compared to those in solution (8). By comparing all six types of functions involved in solutions (7) and (8), the functions can be ranked by their degree of nonlinearity, from the weakest to strongest: linear, logarithmic, inverse time, exponential, power-law and exponential with inverse time.

When analyzing both the singularity and nonlinearity of the six functions as t → 0, we can combine their behaviors into a comprehensive framework. The details of nonlinearity and singularity of these functions are summarized in

Table 1. Nonlinearity and singularity analysis of functions as $\mathrm{t}\to 0$.

Function	Form of the Function	Singularity	Nonlinearity	Combined Behavior
Linear Function	f(t) = t	Non-singular	None	Smooth, finite and entirely predictable.
Logarithmic Function	f(t) = log(t)	Weakly singular	Weak	Exhibits mild nonlinearity with weak singularity.
Inverse Time Function	f(t) = 1/t	Strongly singular	Moderate	A strongly singular and moderately nonlinear function.
Exponential Function	f(t) = exp(−λ0t)	Non-singular	Moderate	A non-singular function with moderate nonlinearity.
Power-law Function	$\mathrm{f}\left(\mathrm{t}\right)=t^{-{\lambda }_1}$	Moderately singular	Strong	A strongly nonlinear and moderately singular function.
Exponential with Inverse Time Function	f(t) = exp(λ2/t)	Strong singularity	Very strong	A strong singular and highly nonlinear function.

.

Linear functions are the most straightforward, with no singularity or nonlinearity. Logarithmic functions diverge weakly and are only slightly nonlinear. Inverse time functions (special case of power-law function with −1 as exponent) exhibit both strong singularity and moderate nonlinearity, making them more sensitive near t = 0. Exponential functions are non-singular but moderately nonlinear, offering a smooth and finite decay. Power-law functions are strongly nonlinear and moderately singular, with divergence strength controlled by λ₁. Exponential with inverse time functions are strongly singular and exhibit strong nonlinearity, combining rapid decay with sensitivity near t = 0. Figure 1 illustrates the curves of the linear, logarithmic, exponential and power-law functions fitted to a same dataset for visualization purpose.

This combined analysis reveals the interplay between singularity and nonlinearity, with some functions (e.g., power-law) showing both strong properties, while others (e.g., exponential) balance nonlinearity with smoothness. Model (6) can be used as a combined model as used in [12] for describing the frequency of floods in the Great Toronto Area, in [13] for describing the frequency distribution of the grade of porphyry copper deposits of the world, and in [14] for describing geochemical anomalies in a vertical profile in a mineral district in Inner Mongolia. The combined model of power-law and power-law with exponential cutoff was also applied in [4] for biodiversity estimation and turning-point prediction for Coronavirus Diease-2019 (COVID-19) infection/fatality. One property that distinguishes the power-law model from other functions is its scale-invariance or self-similarity, which is commonly associated with fractal systems [15,16]. To model irregularities around a location or time point with a singularity in density, we turn to fractal or fractional calculus, where derivatives are defined at non-integer orders. For further discussion on fractal density and fractal calculus, refer to [17].

3. Magmatic Activities Concerning Tectonic Evolution

The NCC, one of the oldest cratons on the globe, is composed of the Eastern Block, the Western Block and the Trans-North China Orogen (TNCO) [18]. The Western Block comprises the Yinshan Block, Orodos Block and the Khondalite Belt. The TNCO is further subdivided into seven metamorphic terranes from north to south: Chengde, North Hebei, Xuanhua, Huai’an, Hengshan, Wutai, Fuping, Lüliang, Zanhuang, Zhongtiao, Dengfeng and Taihua terranes [18] (Figure 2). The study area selected for this research is the Huai’an terrane, situated on the northwestern edge of the TNCO, at the conjunction of the TNCO and the Khondalite Belt of the Western Block in the NCC (Figure 2). The Huai’an terrane has undergone complex tectonic evolution and intense magmatic activity during the Paleoproterozoic, making it an important window for understanding the tectonic processes of the NCC [19,20,21].

The tectonic evolution of the NCC during the Paleoproterozoic can be divided into two primary phases of activity. In the first phase, the amalgamation of seven micro-blocks around 2.5–2.3 Ga formed a stable cratonic block, marking the initial cratonization of the NCC [22,23]. The second phase of cratonization occurred around 2.3–1.82 Ga, involving continental rifting (2.35–1.95 Ga) and subduction–accretion–collision processes (1.95–1.82 Ga) [24,25,26,27]. During 1.95–1.80 Ga, the Huai’an terrane underwent three episodes of metamorphism and deformation, which are represented by crustal thickening, extrusion and exhumation, respectively [28]. These three episodes of metamorphism and deformation in the Huai’an terrane during 1.95–1.80 Ga indicate the NW-SE compressional tectonic regime, marking the final formation of the NCC [29,30]. Previous research about the Huai’an terrane has primarily focused on the tectonic background during the Late Archaean to Early Paleoproterozoic, as well as the metamorphic events and magmatic events during the Late Paleoproterozoic [31,32,33,34,35,36]. However, there is a lack of in-depth studies on the singularities of magmatic events concerning crustal growth processes during the Middle to Late Paleoproterozoic [37,38].

3.1. Datasets

To understand the magmatic events and tectonic evolution of the Huai’an terrane, 12 granite samples were collected along the Datong–Shangyi fault zone, with nearly equal spacing from north to south of the Huai’an terrane for zircon U-Pb dating (Figure 2). The 129 spots analyzed for zircon U-Pb ages of these 12 granites are listed in Supplementary Table S1. Additionally, 73 granite samples in the TNCO from 26 published papers were incorporated to provide a more comprehensive dataset for zircon U-Pb age analysis. The 1439 spots analyzed for zircon U-Pb ages of these 73 granites are listed in Supplementary Table S2. The combined datasets are used herein to investigate magmatic activity and perform singularity analysis of age peaks, offering valuable insights into the timing and nature of magmatic events and their relationship to tectonic processes in the region.

The combined dataset includes zircon age data from various magmatic events in the TNCO, and their statistical parameters are shown in

Table 2. Statistical parameters of granites’ magmatic zircon ages from the three phases of magmatic activity in the TNCO.

Age Range (Ga)	2.6–2.2	2.2–1.9	1.9–1.6
Number of Samples	528	530	484
Mean (Ma)	2448.30	2076.51	1812.73
Median (Ma)	2473	2075	1820
Standard Deviation	95.42	68.43	48.43
Skewness	−0.14	−0.12	−0.91
Kurtosis	0.12	0.14	1.25

. This dataset comprises zircon ages from three distinct magmatic phases, with the age ranges 2.6–2.2 Ga, 2.2–1.9 Ga and 1.9–1.6 Ga. The statistical analysis includes the number of samples, mean age, median age, standard deviation, skewness and kurtosis for each of the three age ranges. These statistics provide insight into the distribution and variability of zircon ages, helping to characterize the magmatic activity in the region.

The granites’ zircon magmatic age data are presented in the histogram in Figure 3. The ages of granites in the TNCO mainly fall into three intervals, 2.6–2.2 Ga, 2.2–1.9 Ga, and 1.9–1.7 Ga, with peak values occurring around ~2.53 Ga, ~2.08 Ga, and ~1.84 Ga, respectively. These age distributions highlight key magmatic phases in the region, and the spatial maps in Figure 2 illustrate the geographical extent of these magmatic events.

3.2. Analysis of Magmatic Activity Using Nonlinear Functions

To analyze the variations in magmatic activity, several functions derived in this study in Section 2 were employed, including power-law, logarithmic and exponential functions, each chosen to capture different aspects of nonlinearity and singularity in the data. These functions were applied to quantify the behavior of magmatic activity around the three peaks of the magmatic events.

Statistical analysis was performed using age bins of 20 Ma, 30 Ma and 40 Ma intervals, and the data were converted into time series to represent the probability distribution, minimizing biases due to inconsistent bin sizes. For each time series, a density of the number of age data per unit tie interval around a peak was calculated, and the nonlinear functions derived in this study were applied for analysis. The age density of the time series is calculated as follows:

$$\rho (\Delta t)={\displaystyle \frac{1}{\Delta t}}N(\Delta t),$$

(9)

where $\Delta t$ represents the age interval around an age peak, $N(\Delta t)$ is the number of zircon ages within the age range of t ± 1/2Δt and $\rho (\Delta t)$ represents the average age density [16].

The average age densities for both sides and one side of the three magmatic event peaks were calculated. One of the six types of functions as shown in Table 1 was fitted using the least squares method. The left side of the peak corresponds to younger ages, while the right side corresponds to older ages. The analysis results for the three peaks are displayed in Figure 4, Figure 5 and Figure 6.

The results indicate that the trends on both sides of the peaks for the three magmatic episodes are similar. However, the density attenuation curves fitted by the LS method exhibit significant differences in form, particularly showing evident asymmetry around the peaks.

For the magmatic activity during 2.6–2.2 Ga, the average attenuation variations on both sides of the peak (2.53 Ga) can be fitted with logarithmic functions, indicating a logarithmic decrease in anomaly intensity from the center outward. On the right side of the peak (earlier than 2.53 Ga), the density of magmatic activity increases rapidly, following a power-law form, while after the peak, the activity decreases logarithmically. This asymmetry suggests that the magma intrusion rate was fast, while the decline and cessation of magmatic activity were relatively slow. The earlier period (log-log linear) exhibits stronger nonlinearity compared to the later period (log-linear).

For the magmatic activity during 2.2–1.9 Ga, the magmatic activity on both sides of the peak (earlier than 2.08 Ga and later) follows either logarithmic or exponential attenuation functions, indicating a less pronounced asymmetry compared to the earlier period.

For the magmatic activity during 1.9–1.6 Ga, the attenuation around the peak (1.84 Ga) and on both sides is characterized by exponential functions with small exponents, reflecting a weaker degree of nonlinearity compared to the earlier two episodes. This final episode no longer exhibits singularity.

The three magmatic episodes exhibit distinct nonlinear and singularity characteristics, with strong linearity and singularity represented by a power-law function fitted to the first peak on the right side, transitioning to a mix of logarithmic and exponential functions in the second peak, and culminating in purely exponential functions at the last peak. These variations likely reflect differences in the deep mechanisms driving the magmatic activity as to be discussed.

4. Discussion

4.1. Magmatic Events and Crustal Growth and Reworking

Geochemical and isotopic data of granitoid rocks from the Chengde and Wutai terranes highlight the widespread magmatic activity during 2.6–2.2 Ga [39,40,41]. Granitic rocks of 2.6–2.54 Ga indicate crustal growth and magmatic processes associated with the assembly of the early continental crust. Specifically, early Proterozoic granitic magmatism is well documented in the TNCO, with zircon εHf(t) values and T_DM data suggesting a relatively uniform magma source for 2.2–1.9 Ga granitoids [42]. The presence of A-type and I-type granites in Dengfeng and Hengshan terranes (2.50 Ga and 2.37–2.47 Ga, respectively) points to crustal growth through mantle-derived magmatism. The TNCO during 2.6–2.2 Ga is in a subduction tectonic setting, where the subduction of the oceanic crust beneath the continental crust led to mantle wedge melting and the intrusion of basaltic magma into the lower crust, further increasing the overall volume of newly generated crust [38,43,44] (Figure 7a). This mantle-driven magmatism is consistent with the evolution of an active continental margin, forming the foundation for the crustal growth of the TNCO during 2.6–2.2 Ga.

The geochemical data from various rocks in the TNCO indicate a complex transition from crustal growth to reworking during 2.2–1.9 Ga [35]. In the Hengshan terrane, the 2.05–2.06 Ga granite exhibits negative εHf(t) values, suggesting that the magma primarily derived from the partial melting of the lower crust, indicating crustal reworking [45]. The 2.14 Ga granite from the Wutai terrane displays mixed isotopic signatures, indicating a combination of mantle and lower crustal sources, signifying crustal growth and reworking [46]. In the Fuping terrane, the A-type granites at 2.08–2.07 Ga show positive εHf(t) values, indicating that the magma originated from newly formed, depleted mantle or lower crust, which reflects crustal growth and reworking [47]. These findings suggest that during 2.2–1.9 Ga, the TNCO transitioned from crustal growth to reworking, with the involvement of mantle-derived and ancient lower crustal materials. This transitional phase was influenced by tectonic activities such as subduction and continental rifting, with notable mantle-derived basaltic and acidic magmatic events marking this era [48,49].

Geological evidence from the terranes of Chengde, Lüliang and Taihua indicates extensive granite magmatism during 1.9–1.7 Ga, pointing to the partial melting of ancient lower crust as the source of these magmas [38,39,50]. The εHf (t) values of zircon grains from the 1.9–1.7 Ga granite magmatism show considerable variability, suggesting complex sources and heterogeneous composition of the underlying crust [51]. Specifically, in the Chengde terrane, the 1.91–1.75 Ga granitic magmatism indicates the reworking of the ancient lower crust, with the εHf (t) values ranging from −6.20 to −1.70 and T_DM1 values between 2.33 and 2.50 Ga, supporting the idea of crustal remelting [52]. The Lüliang terrane also exhibits evidence of crustal reworking through partial melting of the older lower crust, with magmatic ages ranging from 1.88 to 1.79 Ga [50,53]. The collision between the Eastern and Western blocks of the NCC, which likely occurred along the western margin of the TNCO during 1.9–1.7 G, resulted in crustal reworking by the partial melting of older lower crust (Figure 7b).

4.2. Nonlinearity and Singularity of Magmatic Events

The differences in nonlinearity and singularity across the three magmatic periods (2.6–2.2 Ga, 2.2–1.9 Ga and 1.9–1.6 Ga) can be linked to the evolving tectonic settings and processes driving magmatic activity during each period:

(1)

Period 1: 2.6–2.2 Ga (Late Archaean to Early Proterozoic)

The power-law behavior of the density function formed prior to the peak (rapid increase in magmatic activity) suggests a fast, dynamic process, such as the rapid intrusion of magma due to active subduction zones. The transition to logarithmic decay after the peak reflects a gradual cessation of activity, as magmatic processes slowed down. This strong nonlinearity is indicative of a tectonically driven burst of magmatic activity, driven by the rapid convergence and subduction processes. The singularity on the right side of the peak could indicate a highly focused, intense magmatic episode, such as crustal growth via mantle-derived magmatism, with sudden onset and a slow tapering off. During 2.6–2.2 Ga, the strong nonlinearity (power-law function and logarithmic function) and singularity are linked to the dynamic and intense subduction-related magmatism, where rapid mantle-derived intrusion and crustal growth led to quick magma emplacement and slower cessation.

(2)

Period 2: 2.2–1.9 Ga (Early Proterozoic)

The less pronounced asymmetry and the shift to logarithmic or exponential attenuation functions indicate a less intense but more stable magmatic environment. The decrease in nonlinearity reflects a tectonic environment where magmatic activity was still driven by subduction or rifting but with a more gradual and less explosive nature. This indicates that tectonic forces had stabilized, and the magma generation was less sudden, with more gradual and prolonged magmatic events. During 2.2–1.9 Ga, the less pronounced asymmetry and more gradual decay (logarithmic/exponential) reflect the transition phase from crustal growth to reworking and complex tectonic interaction (e.g., rifting and arc-continent collision) resulted in less intense but prolonged magmatic activity.

(3)

Period 3: 1.9–1.6 Ga (Paleoproterozoic)

The weaker nonlinearity and the shift to purely exponential decay after the peak suggest a much more mature and stabilized tectonic setting. Magmatic processes during this period were likely more uniform and continuous, with fewer sudden peaks in activity. The lack of singularity indicates a more homogeneous magmatic system and the absence of large, sudden influxes of magma from the mantle. Instead, magma derived from crustal partial melting led to a more gradual, steady decay in magmatic activity. During 1.9–1.6 Ga, the weak nonlinearity and absence of singularity indicate a period of stable tectonics with crustal partial melting. The more uniform and continuous magmatic processes reflect a tectonic environment where continental crusts had stabilized and magmatic activity was much less intense.

In conclusion, these differences between the three magmatic events are tied to the changes in tectonic regimes, where the earlier periods were more dynamic with rapid mantle processes, while later periods reflected the slow, steady evolution of crustal structures.

5. Conclusions

This study highlights the application of nonlinear models—specifically power-law, logarithmic and exponential functions—in understanding the temporal and spatial dynamics of magmatic events, with a particular focus on the Huai’an terrane. Through the analysis of U-Pb zircon dating, Lu-Hf isotopes and geochemical data, the research has shown that these nonlinear models offer valuable insights into the complex tectonic evolution of the NCC, revealing distinct phases of crustal growth and reworking. The power-law function and logarithmic functions effectively capture the intense magmatic activity of 2.6–2.2 Ga, with clear singularity and high nonlinearity. In contrast, the mix of logarithmic and exponential models with small exponents describe the transitional phases (2.2–1.9 Ga and 1.9–1.7 Ga), where magmatic processes evolved more gradually due to changing tectonic conditions, such as subduction and continental rifting.

By applying the nonlinear models, this study provides a new framework for analyzing extreme geological events, emphasizing the nonlinear nature of geological processes. This work may contribute to the broader field of tectonic evolution by demonstrating how nonlinear functions, particularly those exhibiting singularity behaviors, can be used to unravel the intricate dynamics of Earth’s magmatic history. It must be remembered that the degree of statistical significance of LS fittings to the data with different functions is relatively small. More high-quality data are needed to validate the findings in the study area or other areas in further study. Given the quality of data and the limited case studies focused on a single area within the NCC, the results may not be broadly generalized to other regions or events. Further studies and validations are required to strengthen this line of research. Nevertheless, the distinctions and similarities among the nonlinear models discussed, as well as their quantification based on nonlinearity and singularity, offer valuable insights. These findings can facilitate the effective utilization and interpretation of these functions in describing anomalies associated with extreme geological events.