Modeling forest soil organic carbon dynamics under seasonal freeze–thaw cycles: implications for carbon sequestration in cold-region forest ecosystems

In temperate to cold-temperate forest regions, seasonal freeze–thaw processes exert a significant influence on soil organic carbon (SOC) dynamics and carbon balance. These regions act as critical carbon sinks, where low temperatures and periodic freezing inhibit organic matter decomposition. However, under global warming, the shortening of frozen periods and the intensification of freeze–thaw cycles may accelerate the decomposition of soil organic matter, releasing more carbon into the atmosphere. Accurately characterizing SOC variations under freeze–thaw conditions is therefore essential for understanding terrestrial carbon cycling and improving regional carbon assessments. In this study, a typical seasonally frozen forest area in the Xiaoxing’an Mountains of northeastern China was selected as the research site. A hybrid deep learning framework integrating Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN), Long Short-Term Memory (LSTM) networks, and two-dimensional Convolutional Neural Networks (CNN) was developed to simulate the temporal dynamics of forest SOC. The CEEMDAN method was used to decompose the non-stationary SOC sequence into multiple intrinsic components, effectively improving the model’s ability to capture nonlinear and multi-scale variations associated with freeze–thaw cycles. The CEEMDAN–LSTM–CNN model demonstrated substantially higher predictive accuracy than conventional CNN–LSTM and LSTM models, with R 2 increased by 37.6 and 47.8%, respectively, and error metrics (RMSE, MAE, MAPE) reduced by more than 30%. The optimized model achieved a stable R 2 of 0.83, indicating strong predictive performance. These results show that the proposed hybrid framework provides a robust and accurate approach for modeling SOC dynamics in seasonally frozen forest soils, supporting improved understanding of carbon processes and their response to climate change in cold-region ecosystems.

Permafrost and seasonally frozen ground regions are long regulated by cold and humid environmental conditions, where the decomposition rate of soil organic matter is markedly suppressed, leading to extensive accumulation of soil organic carbon (SOC). These regions thus constitute one of the most critical carbon reservoirs in the global terrestrial ecosystem (Plaza et al., 2019; Zimov et al., 2006). Existing studies estimate that soils in the Northern Hemisphere’s permafrost and seasonally frozen regions store more than 50% of the world’s total SOC, with a carbon stock approximately twice that of the current atmospheric carbon pool (Zimov et al., 2006).

The major carbon reservoirs in these frozen regions highly overlap with the forest ecosystems of the mid- to high-latitude Northern Hemisphere, tightly linking the stability of permafrost carbon pools with the carbon sink functions of forests. However, intensifying global warming and increasing human disturbances are significantly undermining the stability of these frozen environments (Kong et al., 2023; Wu et al., 2021; Rogelj et al., 2012). The growing sensitivity of permafrost systems to external perturbations not only poses a severe threat to the long-term sequestration capacity of these vast carbon pools but also triggers cascading effects on the sustainability of forest carbon sinks in overlapping regions (Jorgenson et al., 2001; Yang et al., 2019).

Under these conditions, forest ecosystems in the mid- to high-latitude Northern Hemisphere, which are key components of the terrestrial biosphere’s carbon sinks, have become a focal frontier in global change research. This is particularly true for the sustainability of their carbon sequestration functions (Bonan, 2008; Nave et al., 2024; Wang et al., 2022). Multi-scale observational data indicate that, driven by both climate warming and human activities, the carbon cycling processes of forest soils in seasonally frozen regions are undergoing significant shifts, manifested by declining soil fertility and continuous reduction of surface SOC. This trend signals a potential weakening of the regional carbon sink capacity (Zhang et al., 2021; Gao et al., 2021).

However, the underlying drivers of these observed trends remain incompletely understood, and elucidating these mechanisms is essential for predicting future changes in carbon sink capacity. Mechanistically, the dynamic evolution of soil carbon pools is governed by complex feedback interactions between the climate and soil systems (Crowther et al., 2019). If warming-induced SOC mineralization continues to exceed the carbon input from vegetation net primary productivity, forest ecosystems may shift from acting as carbon sinks to becoming carbon sources. This transformation not only disrupts carbon balance at regional and global scales but also intensifies global warming through positive carbon and climate feedback loops, leading to a “warming–carbon release–further warming” cycle (Davidson and Janssens, 2006; Schimel et al., 1994; Shao et al., 2023).

Notably, freeze–thaw cycling, which is a core physical process in high-latitude and high-altitude regions, is highly sensitive to temperature variations (Grogan et al., 2004; Li et al., 2024; Yin et al., 2022; Mayr et al., 2020). Under warming scenarios, alterations in freeze–thaw rhythms are expected to weaken the physical protection of organic carbon by low temperatures, accelerating the decomposition of historically sequestered carbon, thereby reducing soil carbon storage and amplifying the positive feedback of the carbon cycle to climate change (Kirschbaum, 1995; Wang et al., 2025).

Consequently, elucidating the dynamic response mechanisms of forest SOC under seasonal freeze–thaw effects and accurately quantifying the threshold between carbon sink and carbon source transitions have become core scientific challenges for simulating future global carbon cycle trajectories, assessing climate feedback intensities, and formulating effective climate mitigation strategies.

However, most existing studies investigating SOC dynamics in mid- to high-latitude freeze–thaw regions have relied on traditional multiple linear regression methods, which are limited in capturing complex nonlinear relationships (Zhu et al., 2025; Bortolo et al., 2025; Shi et al., 2025). These methods generally assume linear and stationary temporal dependencies, limiting their ability to characterize the nonlinear fluctuations and nonstationary patterns typical of SOC dynamics under freeze–thaw conditions. For example, Hengl et al. (2015) compared random forest and linear regression approaches for predicting soil properties, including soil organic carbon, across the African continent. Their results showed that random forest models consistently achieved higher prediction accuracy than linear regression models, with lower prediction errors across different soil properties and depths, highlighting the limitations of linear models in representing complex soil property variations. Therefore, for SOC dynamics under seasonal freeze–thaw conditions, where temporal variations often exhibit nonlinear fluctuations and nonstationary characteristics, more flexible nonlinear modeling approaches are required to improve prediction accuracy. In recent years, deep learning models such as Long Short-Term Memory (LSTM) and Convolutional Neural Networks (CNN) have been increasingly applied to spatiotemporal data modeling and have demonstrated significant advantages in capturing complex nonlinear patterns (LeCun et al., 2015; Zhang et al., 2022).

The CNNs perform well in processing structured data by optimizing memory use, reducing parameters, and mitigating overfitting. Their ability to extract local features can also be extended to time-series analysis, enabling the identification of critical patterns in SOC data to improve predictive accuracy (Padarian et al., 2019). However, SOC dynamics in seasonally frozen regions exhibit strong non-stationarity, abrupt transitions, and multi-scale fluctuations caused by freeze–thaw cycles, climate variability, and coupled soil–vegetation processes. Although CNNs can effectively capture local temporal features, their convolution-based structure may have limited capability in representing long-range nonlinear interactions and evolving temporal dependencies associated with freeze–thaw-driven SOC variations.

LSTM networks, through their gating mechanisms, effectively overcome the vanishing and exploding gradient problems inherent in traditional Recurrent Neural Networks (RNN), allowing them to capture long-term dependencies and utilize historical information to enhance present predictions. They have become essential tools for nonlinear modeling of high-dimensional and complex time series (Wang et al., 2022). Nevertheless, LSTM models may be sensitive to high-frequency noise and irregular fluctuations commonly present in SOC observations, which can originate from field sampling uncertainty, spatial heterogeneity of soil properties, and environmental disturbances. These factors may obscure meaningful temporal signals and reduce model stability and interpretability under freeze–thaw conditions.

Despite CNN’s capacity to process spatially structured inputs via convolution, pooling, padding, and batch normalization, it remains limited when addressing the highly nonstationary SOC dynamics induced by freeze–thaw processes. Conversely, traditional temporal models such as LSTM are sensitive to high-frequency noise, which may obscure critical signals and diminish the model’s interpretability under freeze–thaw-driven mechanisms.

To address these issues, Huang et al. innovatively proposed the Empirical Mode Decomposition (EMD) method, which decomposes time-series signals into a series of Intrinsic Mode Functions (IMFs) representing different scales of trend or fluctuation, thereby enabling the decomposition of nonstationary and nonlinear signals (Huang et al., 1998). To mitigate the problems of incomplete decomposition, false components, and mode mixing in EMD, Wu et al. introduced uniformly distributed auxiliary noise to develop the Ensemble Empirical Mode Decomposition (EEMD) method, which resolved mode mixing but left residual Gaussian white noise (Wu and Huang, 2004). Torres et al. further refined EEMD by adding adaptive white noise, proposing the CEEMDAN method.

CEEMDAN effectively eliminates incomplete decomposition and large reconstruction errors, producing relatively simple and independent IMFs that facilitate the extraction of key fluctuation features (Torres et al., 2011). Compared with traditional signal decomposition methods, such as Fourier transform and wavelet analysis, CEEMDAN does not require predefined basis functions and can adaptively extract intrinsic oscillatory components from nonlinear and nonstationary signals. Compared with EMD and EEMD, CEEMDAN further reduces mode mixing and residual noise through adaptive noise-assisted decomposition, thereby improving signal reconstruction quality and prediction reliability.

Previous studies have demonstrated the effectiveness of CEEMDAN in environmental time-series forecasting, including applications in climate variables, hydrological processes, and hydro-climatic analysis (Sayadi et al., 2025; Jayaparvathy et al., 2024; Li et al., 2024). These successful applications indicate that CEEMDAN provides an effective approach for extracting complex temporal patterns from environmental datasets characterized by nonlinear fluctuations and multiple frequency components.

Existing soil time-series prediction studies primarily use LSTM and its extended architectures (e.g., CNN-LSTM) as core modeling frameworks (Azad and He, 2025; He et al., 2025; Farhangmehr et al., 2025). However, these models generally assume data stationarity, which is inconsistent with the markedly nonstationary nature of SOC dynamics in freeze–thaw regions influenced by multiple coupled processes. This inherent mismatch constrains their capacity to describe SOC evolution under freeze–thaw effects. Additionally, observed SOC data are often compromised by various disturbances such as field sampling bias, spatial heterogeneity in soil physicochemical properties, and instrumental errors, leading to unavoidable noise and random fluctuations that degrade both prediction accuracy and model stability.

To address these methodological and data-related issues, a hybrid CEEMDAN-LSTM-CNN predictive model was constructed in this study by integrating CEEMDAN-based signal denoising with an LSTM-CNN time-series modeling framework. The CEEMDAN component adaptively decomposes SOC sequences into trend, periodic, and noise components, while the LSTM-CNN structure leverages LSTM’s temporal dependency learning and CNN’s local feature extraction capabilities. Together, they form a synergistic framework that enhances SOC prediction performance through improvements in both data preprocessing and model construction.

The specific objectives of this study are threefold: (1) to develop a hybrid deep learning framework for improving the prediction accuracy of forest SOC dynamics under seasonal freeze–thaw conditions; (2) to enhance the extraction of nonlinear temporal features by integrating CEEMDAN-based signal decomposition with LSTM-CNN architectures; and (3) to provide a reliable methodological framework for understanding SOC responses and supporting forest carbon management in climate-sensitive regions. By combining adaptive signal decomposition and deep learning approaches, this study aims to overcome the limitations of conventional statistical and single-model approaches and provide new insights into forest carbon-cycle responses under future climate change scenarios.

Study area and data

Overview of the study area

This study focuses on a natural mixed coniferous–broadleaf forest in the Xiaoxing’an Mountains, northeastern China. The experimental site was located in the Dongfanghong Forest Farm, Dailing Forestry Bureau, Heilongjiang Province. The forest farm is located approximately 13.9 km southeast of Dailing District, within the coordinates 46°50′8″–47°21′32″N and 128°37′46″–129°17′50″E (Figure 1).

The region represents a typical seasonally frozen area, where the soil begins to freeze in late October and remains frozen until mid-to-late April of the following year, showing pronounced freeze–thaw cycling. The mean annual temperature is 2.7 °C. Winters are cold, dry, and snowy with weak winds, while summers are cool, humid, and rainy. The site has an average slope of 10°, a mean elevation of 504 m, and an average annual precipitation of 661 mm.

Soils are primarily classified as dark brown forest soils, with small patches of meadow and marsh soils found in valley areas. The average soil thickness is approximately 40 cm. The forest community at the sampling site is a natural mixed coniferous–broadleaf forest. The dominant tree species include Korean pine (Pinus koraiensis), fir (Abies fabri), spruce (Picea asperata), linden (Tilia tuan), mono maple (Acer mono Maxim), David’s maple (Acer davidii Franch), and Manchurian ash (Fraxinus mandshurica Rupr). The average stand age is about 70 years, with a canopy closure degree exceeding 0.8.

Experimental design and soil sampling

In this study, four representative periods were selected to conduct field sampling and to characterize the annual freeze–thaw dynamics of seasonally frozen forest soils, namely the unfrozen period (July), the freeze–thaw period (October), the frozen period (January), and the thawing period (March) (Figure 2). Within the designated standard site, seven sampling plots were established. Each plot contained three survey quadrats of 30 m × 30 m, and within each quadrat, two 1 m × 1 m soil sampling units were randomly arranged. Before sampling, the surface litter layer was carefully removed. Using a 20 cm × 20 cm sampling frame, surface soil samples were collected from a depth of 0–10 cm. This soil layer represents the most biologically active part of the forest soil carbon cycle and is most sensitive to freeze–thaw processes (Kuzyakov et al., 2000; Mellander et al., 2007). After collection, all soil samples were immediately transported to the laboratory to serve as the basic materials for subsequent temperature-controlled incubation experiments.

The temperature settings for the incubation experiments were determined according to field-monitored annual meteorological data. The 15-day mean air temperatures corresponding to the four representative periods were used as benchmarks: 17 °C–30 °C for the unfrozen period, −5 °C–5 °C for both the freeze–thaw and thawing periods, and −30 °C to −16 °C for the frozen period. All treatments followed an alternating temperature cycle, with high and low temperatures switched every 12 h to form a complete 24-h cycle. A total of 30 cycles were conducted to simulate repeated freeze–thaw events.

Soil samples were collected at the end of the 0th, 1st, 3rd, 5th, 7th, 15th, and 30th cycles. To avoid repeated sampling interference and ensure the accuracy of soil property measurements, independent incubation boxes were used for each sampling time point. In total, 1,176 incubation boxes were used (4 periods × 7 plots × 2 sampling units × 7 time points × 3 replicates). After incubation, soil samples were air-dried, ground, and sieved through a 0.15 mm mesh. The SOC content was determined using the external heating potassium dichromate–sulfuric acid oxidation method, while the initial SOC content was directly measured from the field-fresh samples.

Modeling algorithm

Long short–term memory network (LSTM)

The LSTM network is a specialized form of RNN developed to capture long-range dependencies in sequential data (Hochreiter and Schmidhuber, 1997). It effectively addresses the vanishing or exploding gradient problems common in traditional recurrent neural networks (RNNs) during training. This is achieved through its unique gating mechanisms and a directed update strategy for the cell state. The core architecture of an LSTM unit includes a Forget Gate, an Input Gate, an Output Gate, and a Cell State for retaining long–term information (Figure 3). The mathematical principle of the LSTM model is formulated in Equations 1–6:

The function of the forget gate is to regulate the extent to which the prior cell state is preserved. This is achieved through the output f t, which is generated by processing the current input x t and the previous hidden state h t − 1 with their corresponding weight matrices U f and W f, and adding the bias term b f.

The input gate regulates the extent to which the current input information updates the cell state. i t is the activation output of the input gate, and C t ˜ denotes the candidate cell state. W i 、 U i 、 W c and U c are the weight coefficients for the linear relationships, while b i and b c are bias terms.

Subsequently, the cell state is updated. C t represents the current cell state, and C t − 1 is the cell state from the previous time step. The symbol ⊙ indicates element–wise multiplication.

The output gate generates the hidden state based on the current cell state. O t represents the activation vector, and h t is the current hidden state.

Owing to these gating mechanisms that selectively remember and forget information, LSTMs are adept at capturing long–term dependencies. This capability is crucial for modeling phenomena such as cross–seasonal carbon release dynamics during freeze–thaw cycles. Therefore, in the context of modeling SOC dynamics, the LSTM’s ability to integrate historical states makes it well–suited for capturing temporal dependencies inherent in complex ecological processes.

Convolutional neural network (CNN)

Convolutional neural network is a deep learning architecture specialized in processing grid–structured data, such as images, time series, and spatiotemporal sequences. Its key characteristics include local connectivity and weight sharing, which enable efficient extraction of spatial or spatiotemporal features from input data at a relatively low computational cost (Hinton and Salakhutdinov, 2006). A typical CNN structure employs an alternating stack of convolutional layers, pooling layers, and fully connected layers to achieve multi–level feature abstraction from raw inputs:

Convolutional Layer: This layer uses learnable kernels to perform sliding–window operations over local regions of the input data, capturing localized spatial patterns such as edges and textures. The mathematical formulation of the CNN model is given in Equation 7.

Where A i denotes the input, ⊗ represents the convolution operation, σ is the activation function, w i is the weight matrix, and b i is the bias matrix.

Activation Function: A nonlinear function, such as ReLU, is typically applied after the convolution operation. This enhances the model’s ability to represent complex features and helps mitigate the vanishing gradient problem.

Pooling Layer: This layer performs down–sampling operations, such as max pooling or average pooling, to reduce the spatial dimensions of feature maps. This decreases the number of parameters and improves the model’s robustness to transformations such as translation and rotation.

Fully Connected Layer: This layer maps high–level features to the target output space. It connects all nodes from the preceding pooling layer to its own nodes, thereby integrating extracted features and representing the data characteristics effectively.

The core strengths of CNN lie in its local receptive fields, weight sharing mechanism, and hierarchical feature extraction architecture. These properties grant the model translation invariance and scale robustness, making it not only highly effective in image recognition tasks but also widely applicable to local pattern detection in time series analysis (Şener and Tuğal, 2025; Guo et al., 2025).

2.2.3 Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

CEEMDAN is an improved adaptive signal decomposition technique designed to effectively mitigate the mode mixing problem commonly found in traditional Empirical Mode Decomposition (EMD). This method is particularly suitable for multi–scale feature extraction from nonlinear and non–stationary time series. The principle of the CEEMDAN method is presented in Equations 8–11. Its core procedure consists of the following steps:

(1) Adaptive noise injection: A set of adaptive white noise sequences, denoted as δ i (t), is introduced into the original signal X (t) using specific scaling coefficients, generating multiple noise–added signal samples:

Where β₀ is the noise intensity coefficient and N is the number of ensemble trials.

(2) Hierarchical decomposition: Each noise–added signal X i (t) undergoes EMD decomposition. The ensemble average of the first–order intrinsic mode functions (IMF1) is computed as the final IMF1 component:

(3) Residual calculation and iterative decomposition: The first residual is obtained by subtracting the extracted IMF1 component from the original signal. This residual is then treated as the new signal for the subsequent decomposition cycle. The process continues iteratively to derive all subsequent IMF components:

This procedure is repeated until the K–th residual r k (t) becomes a monotonic signal and can no longer be decomposed. The final output includes a series of IMF components IMF ¯ k (t) (k = 1, 2, …, K), arranged in descending order of frequency, along with the final residual r K (t).

Signal reconstruction completeness: The original signal can be exactly reconstructed as the sum of all IMF components and the final residual:

By employing adaptive noise injection and ensemble averaging, CEEMDAN significantly suppresses mode mixing inherent in traditional EMD while ensuring decomposition completeness. This advantage makes it particularly valuable in applications such as time series analysis and multi–scale modeling of SOC dynamics.

Experimental environment configuration

All computations in this study were conducted on a computational platform equipped with a 13th Gen Intel® Core™ i7–13650HX processor (2.60 GHz) and an NVIDIA GeForce RTX 4060 laptop GPU. Neural network modeling was implemented using the TensorFlow–GPU deep learning framework under Python 3.9.22.

The CEEMDAN–LSTM–CNN hybrid prediction model

The CEEMDAN–LSTM–CNN hybrid model proposed in this study is as follows:

(1) CEEMDAN Decomposition.

The CEEMDAN method is employed to adaptively decompose the original SOC time series into a set of IMFs and a residual component, thereby reducing non–stationarity and highlighting the intrinsic temporal modal features within the data.

(2) LSTM–CNN Prediction Module.

A LSTM network is introduced to capture the long–term dependencies and temporal dynamics of each decomposed component. Each IMF component is sequentially input into the LSTM layer (128 memory units, input dimension 10 × 1, return_sequences = True) to learn time–dependent relationships. The LSTM output is then passed through a one–dimensional CNN layer with 512 filters and a kernel size of 3, using the ReLU activation function and ‘same’ padding to extract multi–scale local features. A MaxPooling layer with a window size of 1 is subsequently applied. Dropout regularization (rate = 0.2) is used to prevent overfitting. Finally, the network output is flattened and passed through a fully connected layer to generate the prediction results.

(3) Model Training and Integration.

The model is trained using the Adam optimizer with a learning rate of 0.001 for 100 epochs and a batch size of 100. During prediction, a rolling forecasting mechanism (step size = 10) is adopted. Each IMF component is predicted independently and then linearly recombined, while the residual term is directly inherited as the steady–state component to obtain the final prediction results. The overall model architecture is illustrated in Figure 4.

Comparative models

To validate the effectiveness of the proposed hybrid model, two baseline models were trained for comparison in this study.

(1) LSTM Model.

A time series prediction model consisting of four LSTM layers (Figure 5) was constructed, with input data comprising 15 time steps of historical SOC values. The network structure includes 50, 100, 200, and 300 neurons in successive layers, with the final layer outputting only the state at the last time step. Dropout (rate = 0.2) and a fully connected layer (with 100 ReLU units) are introduced to enhance generalization. The model is trained for 30 epochs using the Adam optimizer with a batch size of 64.

(2) CNN–LSTM Hybrid Model.

This model combines one–dimensional convolution (with 128 filters, kernel size = 1, and ReLU activation function) and bidirectional LSTM (with 64 units) to enhance feature extraction through bidirectional temporal modeling. The model architecture is shown in Figure 6. The network includes multiple regularization techniques such as batch normalization and Dropout (rate = 0.1). The final output is obtained through a fully connected layer. The model is trained for 100 epochs using the Adam optimizer (default learning rate = 0.001) with a batch size of 64.

Model evaluation

To comprehensively evaluate the performance of the prediction models, four widely used evaluation metrics are selected in this study: Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), Mean Absolute Percentage Error (MAPE), and the Coefficient of Determination (R 2). The definitions of these metrics are presented in Equations 12–15:

MAE: This metric reflects the average absolute deviation between the predicted values and the actual values. It is insensitive to outliers, and a lower value indicates higher model accuracy.

RMSE: By amplifying the influence of larger errors through the squared term, RMSE effectively represents the dispersion of prediction errors. A smaller value signifies better model prediction stability.

R 2: This metric measures the proportion of variance in the target variable that is explained by the model. A value closer to 1 indicates a better fit of the model.

Where: m is the number of soil carbon samples; y i ˜ is the actual value of soil carbon at time i; y i is the predicted value at time i; y ¯ is the mean of the actual observed values.

CEEMDAN decomposition and IMF component prediction performance of forest SOC in seasonal freeze–thaw regions

The CEEMDAN–based decomposition method decomposes the original SOC series into 9 IMFs and 1 residue component (Figure 7). From IMF₁ to IMF₉, the decomposition results indicate that the frequency and fluctuation amplitude of the components decrease sequentially. The frequency and amplitude of each component follow a decreasing order, presenting a multi–scale structure from high–frequency details to low–frequency trends. Each IMF component exhibits oscillation patterns that are simpler and more regular than the original sequence, strictly adhering to the theoretical characteristics and convergence conditions of empirical mode decomposition. This confirms the validity of the decomposition results.

Through an in–depth analysis of the physical significance of each IMF component and a comprehensive evaluation of their prediction performance, we can more precisely interpret the underlying driving mechanisms they represent. The prediction performance metrics for each IMF component are shown in Table 1. The high–frequency components (IMF₁–IMF₃) exhibit higher prediction errors, while the low–frequency components (IMF₈–IMF₉) demonstrate significantly lower prediction errors.

Comparison of performance among different prediction models

As shown in Table 2, the proposed CEEMDAN–LSTM–CNN hybrid model demonstrates a significant advantage in the SOC prediction task. Compared to the standalone LSTM model, the RMSE of this model decreased by 32.3%, and the R 2 increased by 47.9%. Compared to the CNN–LSTM model, the RMSE decreased by 22.8%, and the R 2 increased by 37.5%. These results fully demonstrate the effectiveness of combining signal decomposition techniques with deep learning models.

From the perspective of training performance, the accuracy limitations of traditional single models and conventional hybrid models are clearly demonstrated in Figures 8A,B and the corresponding quantitative metrics. As shown in Figure 8A, the fitting curve of the LSTM model exhibits high dispersion, with a determination coefficient (R 2) of only 0.5359. This indicates that when relying solely on temporal dependency for feature learning, the model struggles to effectively capture complex multi–scale fluctuations and key local information within the data, leading to large deviations, high error dispersion, and weak interpretability of the intrinsic data patterns.

In contrast, the CNN–LSTM hybrid model shown in Figure 8B achieves a certain degree of improvement by incorporating CNN’s local feature extraction capability. The RMSE, MAE, and MAPE in the training set are reduced by 2.77, 3.97, and 12.63%, respectively, compared with the LSTM model, while R 2 increases to 0.5621 (a 5.26% improvement). However, due to the lack of multi–scale preprocessing on the raw data, the local features extracted by the CNN still contain redundant noise and irrelevant information, resulting in only limited enhancement and leaving the performance bottleneck of the conventional hybrid architecture unresolved.

The CEEMDAN–CNN–LSTM model, as illustrated in Figure 8C, achieves a qualitative improvement in training performance. Its RMSE, MAE, and MAPE decrease to 7.7199, 5.7507, and 12.411, representing reductions of 36.24, 37.04, and 45.04% compared with the LSTM model, and 34.43, 34.44, and 37.09% compared with the CNN–LSTM model, respectively. The R 2 value rises sharply to 0.8116, improving by 51.45 and 44.39% relative to the LSTM and CNN–LSTM models, respectively. This advantage is reflected in the Taylor diagram in Figure 8G, where the new model’s evaluation index lies closer to the reference point, and in the radar chart in Figure 8H, which highlights its comprehensive superiority across core metrics. The core reason lies in CEEMDAN’s adaptive multi–scale decomposition mechanism, which decomposes the original non–stationary and noisy data into IMFs that capture essential data characteristics. This effectively removes redundant noise and irrelevant fluctuations, providing high–quality, denoised inputs for the subsequent CNN–LSTM network, thereby improving feature fitting efficiency and overall prediction accuracy.

In terms of test set performance, traditional models generally exhibit strong training fits but weak generalization, whereas the new model achieves both high fitting accuracy and strong generalization capability. This distinction is especially evident in the right–hand test set curves of Figures 8A–C and the comparative results shown in Figures 8D–F. For the LSTM model, test set RMSE, MAE, and MAPE increase by 17.45, 19.84, and 16.89% compared with the training set, while R 2 is only 0.5884. Combined with the high dispersion of the test set curve in Figure 8A, this suggests that the model’s reliance on single–sequence dependency results in poor adaptability to unseen noise and data distribution shifts, leading to weak generalization and substantial error fluctuation.

The CNN–LSTM model partially alleviates this issue through CNN’s local feature extraction. Its test set RMSE, MAE, and MAPE decrease by 12.38, 18.23, and 24.12% compared with the LSTM model, and R 2 increases to 0.6702, representing a 13.90% improvement. However, its test set errors remain higher than those of the training set (RMSE +5.84%, MAE + 2.05%, MAPE +1.52%), and the R 2 gap between the training and test sets reaches 10.81%, much larger than that of the new model. The visual discrepancy between the left and right curves in Figure 8B further confirms its limited generalization stability.

By contrast, the CEEMDAN–CNN–LSTM model demonstrates a pronounced advantage in test set performance, with RMSE, MAE, and MAPE of 9.3408, 6.646, and 13.7579, respectively–reductions of 34.32, 39.28, and 47.87% relative to the LSTM model and 25.04, 25.75, and 31.29% relative to the CNN–LSTM model. The R 2 value reaches 0.8186, an increase of 39.12 and 22.14% compared with the LSTM and CNN–LSTM models, respectively. Moreover, the R 2 difference between training and test sets is only 0.7% (compared with 5.25% for LSTM and 10.81% for CNN–LSTM), and the test–train error differences (RMSE 1.6209, MAE 0.8953, MAPE 1.3469) are the smallest among all three, as evidenced by the high overlap between the training and test curves in Figure 8C and the quantitative indices in Figures 8G,H. These results confirm that the new model effectively suppresses overfitting.

In essence, CEEMDAN improves model performance by reducing data non–stationarity and noise at the source through multi–scale decomposition, providing intrinsic feature inputs for modeling. CNN then accurately captures local key features (such as abrupt time–series changes and trend inflection points) within each IMF component, while LSTM explores long–term dependencies across components. This reduces the model’s reliance on noise–driven fitting and enhances generalization. Consequently, the CEEMDAN–CNN–LSTM model significantly outperforms traditional LSTM and CNN–LSTM architectures in both training accuracy and generalization, offering a more reliable modeling framework for complex data–driven prediction tasks.

Synergistic mechanism and generalization capability of the CEEMDAN–LSTM–CNN hybrid model

The core advantage of CEEMDAN decomposition lies in its ability to adaptively decouple complex, non–stationary SOC sequences into a series of IMFs with distinct central frequencies. From the high–frequency IMF₁ to the low–frequency IMF₉ and the residual term, the decomposition results clearly reveal the multi–scale driving mechanisms underlying SOC dynamics. The high–frequency components (IMF₁–IMF₃) capture the most abrupt and frequent short–term fluctuations within the sequence. Their physical origins may stem from random noise, measurement errors, and sudden external disturbances such as short–term heavy rainfall, transient microbial activity bursts, or anthropogenic disturbances. The relatively high prediction errors observed in these components (e.g., IMF₁: RMSE = 9.32, MAPE = 130.22%) do not indicate model failure but rather reflect the inherent randomness and unpredictability of these signals. CEEMDAN effectively isolates these highly uncertain signals from the original data, establishing a foundation for extracting more regular and interpretable information in subsequent analysis.

The medium–frequency components (IMF₄–IMF₇) demonstrate substantially higher predictive accuracy (e.g., IMF₆: MAPE = 10.46%). These components reveal the periodic evolution of SOC at seasonal and interannual scales, likely driven by systematic environmental factors such as temperature and precipitation. Such regular dynamics are particularly well–suited for learning by the LSTM–CNN model. In contrast, the low–frequency and trend components (IMF₈–IMF₉ and Residue) represent the long–term evolution of SOC. Their variations are slow and smooth, reflecting processes operating on multi–year scales–such as soil formation, gradual organic matter accumulation, long–term climate trends, or historical land–use effects. The nearly perfect prediction performance of IMF₈ (RMSE = 0.0099, MAPE = 0.017%) indicates that this component is highly stable and predictable, representing the most deterministic part of the SOC sequence.

In summary, CEEMDAN decomposition serves not only as a signal–processing technique but also as a powerful tool for feature engineering and physical mechanism interpretation. By decoupling a complex original sequence into multiple sub–series at different temporal scales, each corresponding to distinct underlying biogeochemical processes, CEEMDAN transforms the global challenge of modeling SOC dynamics into a set of more targeted sub–problems. This decomposition provides a crucial foundation for achieving high–precision prediction (test set R 2 = 0.814).

The hybrid architecture of LSTM and CNN used in this study achieves complementary strengths in capturing temporal dependencies and extracting local features. The LSTM layer, leveraging its gating mechanism, effectively learns long–term dependencies within each IMF component, capturing both the memory effects and periodic patterns of SOC variations. The subsequent one–dimensional CNN layer functions as a “feature optimizer,” extracting more discriminative local patterns–such as inflection points and abrupt transitions–from the high–level representations generated by the LSTM. This design enhances the model’s sensitivity to subtle dynamic variations, improving both its predictive accuracy and generalization capability.

Optimization of hybrid model architecture parameters and performance regulation

Previous research has confirmed that the CEEMDAN–LSTM–CNN hybrid model achieves significantly higher prediction accuracy than traditional LSTM and baseline CNN–LSTM models. To further enhance model performance, this study focuses on optimizing the LSTM–CNN hybrid architecture through systematic hyperparameter tuning to improve both predictive accuracy and generalization capability. Based on the control variable method, a series of comparative experiments were designed and implemented to investigate the influence mechanisms of each core component on overall model performance. The parameters of the comparative models are listed in Table 3.

The “Base–CEEMDAN–LSTM–CNN” model (No. 1) was used as the benchmark. Its structural configuration includes 128 LSTM hidden units, a single–layer LSTM structure, dropout rates of 0.2 for both the LSTM and CNN layers, 512 convolutional filters, a convolution kernel size of 3, a pooling window size of 1, and a single convolutional layer. On this basis, key network parameters were adjusted sequentially to construct seven variant models, forming eight comparative configurations in total to comprehensively evaluate performance under different structural designs. The learning rate was set to 0.001, with 100 training epochs and a batch size of 100. All experiments were conducted under identical training conditions, data partitions, and evaluation standards to ensure result comparability.

As shown in Figure 9, increasing the number of LSTM hidden units to 256 (No. 2) slightly enhances the model’s nonlinear representation capacity; however, the performance improvement is limited while computational cost increases significantly. This indicates that 128 hidden units already provide sufficient temporal modeling ability for the current task. When the LSTM was further extended to a two–layer stacked structure (No. 3), theoretically it should strengthen deep temporal dependency learning, but in practice it led to higher prediction errors and a reduction in the R 2 value. This suggests that a deeper recurrent structure may cause gradient propagation issues or overfitting, preventing expected performance gains in this task.

In the CNN module, reducing the number of filters to 256 (No.4) weakened the richness of feature extraction, resulting in performance degradation. Similarly, decreasing the convolution kernel size from 3 to 2 (No. 5) limited the receptive field, reducing the model’s ability to effectively identify multi–scale temporal patterns. When the pooling window size was increased from 1 to 2 (No. 6), although spatial invariance and dimensionality reduction efficiency improved, excessive loss of high–frequency detail information occurred, thereby lowering prediction accuracy.

In contrast, introducing a two–layer stacked CNN structure (No. 7) yielded a notable improvement. With other parameters unchanged, constructing a deeper convolutional feature extraction pathway enabled progressive abstraction and efficient fusion of complex local patterns within the input signals. This result confirms the critical role of deepened convolutional structures in feature engineering. Finally, the composite configuration (No. 8)–which simultaneously adopted 256 LSTM units, 256 CNN filters, a kernel size of 2, and a pooling window size of 2–exhibited large performance fluctuations and an overall decline compared with single–parameter optimization schemes. This indicates that simultaneous variation of multiple parameters may disrupt the model’s optimization balance, leading to unstable training or convergence to suboptimal solutions.

Subplots A through D use radar charts to compare the performance of multiple models based on R 2, MAE, MAPE, and RMSE, respectively. Subplot E depicts the influence of hyperparameter tuning on the prediction performance of forest SOC in seasonally frozen regions using a grouped radial bar chart.

In summary, the experimental results reveal that within the CEEMDAN–assisted hybrid modeling framework, the architectural design of the convolutional layer exerts a more pronounced influence on overall model performance. Compared with simply increasing the depth or width of the LSTM, enhancing the hierarchical structure of CNN feature extraction, such as by using stacked convolutional layers, can more effectively improve the model’s expressive power and robustness. Among all configurations, No.7 (LSTM128_StackedCNN_K3_P1) achieved the best performance across all evaluation metrics, confirming the effectiveness and superiority of deepening the convolutional feature extraction pathway. This finding not only highlights the pivotal role of CNN in the hybrid architecture but also provides important empirical evidence and theoretical support for the architectural optimization of deep neural networks targeting complex time series prediction tasks.

Ecological implications, research limitations and future prospects

Although the proposed CEEMDAN–LSTM–CNN framework achieved reliable performance in predicting SOC dynamics under seasonal freeze–thaw conditions, the improvement in prediction accuracy provides important ecological implications. More accurate SOC simulations can enhance assessments of carbon sink stability in cold-region forests, improve understanding of SOC responses to climate warming and freeze–thaw disturbances, and provide scientific support for adaptive forest carbon management and climate mitigation strategies.

Despite these advantages, several limitations should be acknowledged. First, this study was based on soil samples from a representative seasonally frozen forest ecosystem in the Xiaoxing’an Mountains of northeastern China. Therefore, further validation across different forest types, climatic regions, and ecosystem conditions is required to evaluate the generalizability and transferability of the proposed framework. Second, the current model primarily focuses on extracting nonlinear temporal patterns from SOC variations, whereas the underlying ecological mechanisms regulating SOC responses to freeze–thaw processes are not explicitly incorporated. Future studies should integrate key ecological drivers, including microbial activity, soil physicochemical properties, and carbon transformation processes, to improve the mechanistic interpretation of SOC dynamics.

Furthermore, the CEEMDAN–LSTM–CNN framework has potential applicability beyond the studied region. Its adaptive signal decomposition and nonlinear feature-learning capabilities may facilitate SOC prediction in other climate-sensitive ecosystems, including agricultural soils, grasslands, wetlands, alpine ecosystems, and high-latitude forests, where carbon dynamics are strongly influenced by environmental fluctuations and nonlinear processes.

This study proposed a CEEMDAN–LSTM–CNN hybrid model to simulate the dynamic variations of forest soil organic carbon (SOC) in seasonally frozen regions. By integrating CEEMDAN signal decomposition with deep learning, the model effectively addressed the non-stationarity and mode–mixing problems inherent in SOC data. CEEMDAN adaptively decomposed the raw SOC series into intrinsic mode functions, enabling the LSTM and CNN components to capture long–term dependencies and local temporal features, respectively. Experimental results demonstrated that the CEEMDAN–LSTM–CNN model achieved substantially higher prediction accuracy than traditional LSTM and CNN–LSTM models, improving all major performance indicators and enhancing both fitting precision and generalization ability.

Building upon this foundation, further optimization of the CEEMDAN–LSTM–CNN architecture was conducted through systematic hyperparameter tuning. Eight comparative configurations were designed to evaluate the effects of key structural parameters, including LSTM units, convolutional depth, kernel size, and pooling window. Results showed that deepening the convolutional feature extraction path, rather than increasing the LSTM depth or width, most effectively enhanced model robustness and predictive performance. The optimal configuration (No. 7: LSTM128_StackedCNN_K3_P1) achieved the best overall results with a stable R 2 of 0.83, representing a 2.53% improvement over the baseline.

Overall, this study confirms the effectiveness of the decomposition–integration deep learning framework for complex environmental time–series prediction. By combining CEEMDAN–based signal decomposition with hierarchical neural feature extraction, the proposed method provides a robust and transferable modeling approach applicable to hydrological, meteorological, and ecological forecasting tasks. Future research should integrate additional environmental drivers and explore advanced architectures such as Transformer or attention mechanisms to further enhance model generalization and interpretability.