Graph Neural Networks with dilated convolutions for geomagnetic field forecasting: a single-station topological approach

1 Introduction

Geomagnetic field forecasting is a central task in space weather monitoring, as disturbances driven by solar activity can impact technological infrastructure and operational systems. Accurate short-term predictions support mitigation strategies for geomagnetically induced currents, satellite operations, and communication and navigation services. Despite extensive progress in data-driven space weather modeling, improving short-horizon geomagnetic field forecasting remains challenging due to the strongly non-linear and multiscale nature of geospace dynamics Camporeale (2019); Albertson et al. (1974).

The use of machine learning in space weather has grown rapidly, yet forecasting geomagnetic variability remains challenging due to non-stationarity and the multiscale coupling of solar wind drivers with magnetosphere–ionosphere dynamics. Comprehensive reviews emphasize that practical forecasting models must balance predictive skill, interpretability, and strict causality under operational constraints Camporeale (2019). A large fraction of data-driven geomagnetic forecasting has historically focused on global geomagnetic indices (e.g., Dst, SYM-H, Kp), where deep sequence models achieve strong short-horizon accuracy Collado-Villaverde et al. (2021); Abduallah et al. (2020); Collado-Villaverde et al. (2024). However, while global indices provide a macro-scale view of geomagnetic storms, they lack the spatial resolution required for localized risk assessment.

Beyond global indices, forecasting local ground magnetic perturbations is central for applications such as geomagnetically induced current (GIC) risk, where large and rapid magnetic variations are the primary concern. This necessity has motivated community-wide benchmarking efforts, such as the Geospace Environment Modeling (GEM) challenge validation framework, to evaluate model skill under standardized protocols Pulkkinen et al. (2013); Welling et al. (2018). Recent machine-learning studies have revisited these benchmarks with modern deep learning models, actively applying recurrent and convolutional architectures to predict localized ground magnetic perturbations and Geomagnetically Induced Currents (GIC) Pinto (2022). In parallel, large curated magnetometer datasets such as SuperMAG provide standardized preprocessing pipelines that support reproducible station-level modeling and inter-comparison Gjerloev (2012).

Despite these advances, most operational approaches rely on dense magnetometer networks or multi-instrument assimilation. Many practical scenarios, however, require robust single-station forecasts, either due to limited sensor availability, regional monitoring needs, or as building blocks for distributed systems. In this single-station regime, the forecasting problem is deceptively difficult. The three geomagnetic vector components—horizontal ( H ), declination ( D ), and vertical ( Z )—are not independent channels, but coupled manifestations of underlying electrodynamic processes. Although the system is observed at a single geographic location, the underlying state space is inherently structured. This is especially relevant at equatorial latitudes, where ionospheric current systems such as the Equatorial Electrojet (EEJ) strongly modulate the horizontal component and introduce rapid, structured variability Bhattacharyya (2021); Ayebare et al. (2025).

1.1 Spatiotemporal deep learning and research gap

While deep sequence models (like standard LSTMs) have been applied to single-station data, they typically treat the vector components as independent parallel time-series or simply concatenate them, largely ignoring the explicit topological relationship between them. In contrast, Graph Neural Networks (GNNs) provide a principled way to model relational structure beyond Euclidean feature spaces, enabling message passing over a graph to capture directed or undirected dependencies Kipf and Welling (2017); Veličković et al. (2018).

In spatiotemporal forecasting, highly influential architectures combine graph-based propagation with temporal modules to capture both spatial interactions and temporal dynamics simultaneously Li et al. (2018); Yu et al. (2018). A key trend in modern spatiotemporal models is data-driven structure learning: when the true adjacency is unknown or only partially observed, the graph can be inferred from node embeddings. Representative examples include Graph WaveNet and the Multivariate Time-series Graph Neural Network (MTGNN), which explicitly couple graph learning with temporal convolutions for multivariate forecasting Wu et al. (2019), Wu et al.(2020). In parallel, long-range temporal modeling has advanced through causal dilated convolutions and Transformer variants, which achieve large receptive fields while preserving strict causality van den Oord et al. (2016); Bai et al. (2018); Zhou et al. (2021).

However, a significant research gap remains: the application of data-driven structure learning and causal dilated convolutions to the multivariable components of a single isolated station is underexplored. Spatiotemporal models are traditionally designed for physically distributed sensor networks (e.g., traffic grids or multi-station magnetometer arrays). Adapting them to single-station multivariate forecasting requires treating the ( H , D , Z ) components as localized nodes in a latent graph. Furthermore, models operating strictly on single-station ground magnetometer data, without the inclusion of external upstream drivers (such as solar wind parameters or the Interplanetary Magnetic Field, IMF), face an inherent limitation: they cannot fully capture the underlying physical dynamics of the solar-terrestrial coupling. Instead, they must rely on extracting latent topological dependencies and temporal patterns solely from local magnetic perturbations.

To address these gaps, we present a proof-of-concept study for short-term single-station geomagnetic forecasting using 1 year of data (2016) from the MAGDAS network (Jerusalem station, Ecuador) Yumoto and MAGDAS Group (2013); López et al. (2013). By restricting our analysis to this specific spatial and temporal domain, we aim to rigorously evaluate the methodological viability of the proposed architecture. Our objective is to predict the future 10-min evolution of the ( H , D , Z ) vector components by explicitly modeling their cross-component dependencies through a graph-based architecture.

1.2 Contributions

The main contributions of this work are explicitly defined as follows:

• We propose a hybrid Graph Neural Network and Temporal Convolutional Network (GNN-TCN) architecture tailored for single-station forecasting. By modeling the ( H , D , Z ) components as nodes in a learned directed graph, the model infers dynamic spatio-temporal correlations as an adaptive structural regularization, rather than extracting a fixed physical law.
• We provide a comprehensive benchmark against classical Vector Autoregression (VAR) and thoroughly optimized Long Short-Term Memory (LSTM) networks, exploring their respective structural trade-offs. Furthermore, we clarify the operational limits against a naive Persistence baseline, showing that while Persistence benefits from short-term inertia during dominant quiet-time conditions, our model achieves vastly superior phase correlation and tracking fidelity.
• We introduce a spectral characterization of the high-frequency geomagnetic data (using Welch’s Power Spectral Density) to mathematically justify temporal downsampling, demonstrating that fluctuations above the 1-min Nyquist frequency are heavily dominated by instrumental noise.
• We evaluate the model’s robustness not only under continuous data but also stratified by geomagnetic conditions (using a localized dynamic K-index proxy), proving its stability and predictive reliability even during extreme, highly disturbed storm events.

2 Materials and methods

2.1 Data description

The dataset used in this study was obtained from the Magnetic Data Acquisition System (MAGDAS) network, installed by the International Center for Space Weather Science and Education (ICSWSE) at Kyushu University in collaboration with the Quito Astronomical Observatory of the Escuela Politécnica Nacional Yumoto and MAGDAS Group (2013). The station is locally operated by the Quito Astronomical Observatory.

We analyze 1 year of geomagnetic observations (January to December 2016) from the Jerusalem station, located at Jerusalem Park in Malchinguí, Ecuador, which serves as one of the MAGDAS monitoring sites near the geographic equator. Jerusalem is of particular interest due to its proximity to the equatorial region, where ionospheric current systems such as the Equatorial Electrojet (EEJ) strongly modulate the daily variation of the horizontal component Bhattacharyya (2021); Ayebare et al. (2025). The raw MAGDAS station log provides 1 Hz measurements of the three geomagnetic components: Horizontal ( H ), Declination ( D ), and Vertical ( Z ), all measured in nanoTeslas (nT). We model these components as a 3-node multivariate time series.

2.1.1 Spectral analysis and downsampling justification

The raw 1 Hz records were aggregated to a 1-min resolution for multi-horizon forecasting. To mathematically justify this temporal downsampling and characterize the signal-to-noise ratio, we estimated the Power Spectral Density (PSD) of the raw 1 Hz signal using Welch’s method.

The spectral analysis revealed that at frequencies higher than the Nyquist limit corresponding to a 1-min sampling rate (f Nyquist ≈ 8.33 × 1 0 − 3 Hz), the power spectrum flattens significantly. This asymptotic behavior indicates that fluctuations above this frequency threshold are heavily dominated by instrumental white noise rather than underlying physical geomagnetic dynamics. Consequently, aggregating the data to 1-min intervals acts as an essential low-pass filter, effectively removing high-frequency stochastic noise and stabilizing the deep learning optimization process without discarding relevant physical information.

Since data gaps and outages are common in long-term station operations, we applied a cleaning procedure to produce continuous segments suitable for sliding-window learning. The specific resampling, gap-handling, segmentation criteria, and split protocols are detailed in Section 3. Additionally, the data was organized into PyTorch tensors optimized for Compute Unified Device Architecture (CUDA) acceleration. The input structure consists of a temporal sequence of length T i n (representing the historical window of 1-min steps) and 3 spatial channels corresponding to the ( H , D , Z ) magnetic field nodes.

2.2 Methodology

This section details the proposed architecture, which integrates graph structure learning with temporal convolutional networks to model coupled temporal dynamics with latent cross-component structure.

2.2.1 Problem formulation

We frame the monitoring of a single geomagnetic station as a multivariate time series forecasting problem on a graph. Let X ∈ R T × N denote the historical observation matrix, where T is the number of time steps (in minutes) and N = 3 represents the orthogonal components of the geomagnetic field: Horizontal ( H ), Declination ( D ), and Vertical ( Z ).

Although the system is observed at a single geographic location, the underlying state space is inherently structured. Unlike classical statistical models (e.g., Vector Autoregression), which implicitly assume a flat Euclidean feature space where dependencies are inferred solely from global covariance, we model the system as a non-Euclidean domain. Even with N = 3, this formulation explicitly separates which components interact (topology) from how they evolve (dynamics).

We define the underlying topology as a directed, weighted graph G = ( V , E , A ), where:

• V = { v 1 , v 2 , v 3 } is the set of nodes corresponding to the physical state variables { H , D , Z }. Crucially, these nodes represent dimensions in the physical phase space, not spatial coordinates.
• E ⊆ V × V is the set of edges representing functional dependencies. An edge ( v j , v i ) implies that the historical state of component j contains predictive information for the future evolution of component i.
• A ∈ R N × N is the learnable adjacency matrix, where A i j ≥ 0 quantifies the strength of the directional influence j → i. We do not enforce symmetry ( A ≠ A ⊤ ), allowing the model to capture asymmetric causal flows.

Figure 1 illustrates the conceptual shift from the classical flat vector view to our latent graph topology.

At any time step t, the system state is defined as a graph signal x t ∈ R N. Given a historical input sequence of length T i n (in minutes), denoted as X i n = [ x t − T i n + 1 , … , x t ], our goal is to learn a mapping function f ( ⋅ ; θ ) that minimizes the forecasting error over a future horizon of length T out as formulated in Equation 1:where Y ̂ ∈ R T out × N represents the predicted multi-step trajectories for all physical components.

Y ̂ t + 1 : t + T out = f G , X i n ; θ , (1)

2.2.2 Temporal dependency modeling: the TCN module

To capture temporal dynamics without violating the physical causality of the system, we employ a Temporal Convolutional Network (TCN) augmented with Inception-style multi-scale feature extraction.

Unlike Recurrent Neural Networks (RNNs) that process time sequentially, TCNs treat time as a spatial dimension, allowing for massive parallelization while strictly respecting the arrow of time via causal padding.

2.2.3 Dilated convolutions and receptive field

A critical design challenge in geomagnetic forecasting is capturing long-term dependencies without an explosion in parameters. Using standard causal convolutions, the Receptive Field ( R ) grows only linearly with network depth. To address this, we employ causal dilated convolutions, where the filter explicitly skips input steps with a defined spacing factor d. Formally, for a 1D sequence x ∈ R T and a filter w ∈ R K of size K, the operation is defined in Equation 2 as:where d controls the skipping distance between inputs. By stacking layers with exponentially increasing dilation factors (d l = 2 l − 1, i.e., 1,2,4 , …), the receptive field grows exponentially with depth. This mechanism allows the model to integrate information from the distant past with logarithmic complexity. Figure 2 conceptually illustrates how varying the dilation factor d expands the temporal reach of the kernel while strictly preserving causality.

x ∗ d w t = ∑ k = 0 K − 1 w k ⋅ x t − d ⋅ k (2)

2.2.4 Dilated inception layer

Instead of fixing a single kernel size K, which would bias the model towards specific frequencies, our implementation adopts an Inception-style architecture. At each layer, we apply a set of parallel dilated convolutions with distinct kernel sizes K = { 2,3,6,7 } sharing the same dilation factor d. The outputs are concatenated along the channel dimension as shown in Equation 3:this multi-branch design allows the network to simultaneously attend to short-term fluctuations (small K) and broader trends (large K).

H temp = Concat H l ∗ d w k ∣ k ∈ K (3)

2.2.5 Tensor organization

In our pipeline, the input tensor is structured as ( B , C i n , N , T ), where B is the batch size, C i n = 1 is the raw input channel, N = 3 represents the physical nodes (graph vertices), and T the time steps. A 1 × 1 convolution first projects the raw input into a latent feature space ( C out = 32 ). Crucially, the temporal convolutions operate exclusively along the T axis, preserving the independence of the nodes before the subsequent Graph Convolution module mixes them spatially.

2.2.6 Spatial dependency modeling: the GNN module

Classical statistical models often fail to disentangle topology from dynamics, learning a single dense matrix that mixes “who interacts” with “how they interact.” To address this, we propose a Graph Neural Network (GNN) that mathematically decouples structure learning from feature transformation.

The input to this module is a node-feature tensor H ( l ), where the spatial dimension corresponds to the geomagnetic components. We decouple feature transformation from spatial diffusion through a structured three-step process for each layer l:

• a. Temporal Mixing (Dynamics). A temporal convolutional network acts as a feature transformer, mixing information along the time and feature axes independently for each node. This effectively replaces the dense weight matrix of standard GCNs with a temporal inception layer, as expressed in Equation 4:

H temp = ϕ TCN H l (4)

• b. Spatial Mixing (Structure). The learned adjacency matrix A acts strictly as a spatial diffusion operator, propagating the temporally processed features across components:

H spat = A ⋅ H temp (5)

Here, A is a learnable N × N matrix where A i j ≥ 0 (enforced via a ReLU activation function), representing the latent structural influence of component j on i.

• c. Residual Connection. To preserve signal integrity and enable deep stacking without vanishing gradients, we add the input of the block to the output of the spatial mixing step using Batch Normalization (BN), as shown in Equation 6:

H l + 1 = BN H spat + H l (6)

As visualized in Figure 3, these operations enforce a strict mathematical separation: the matrix A learns who interacts, while the shared TCN parameters learn how they interact, ensuring stable temporal feature extraction across components without assuming a static physical law. Furthermore, since the initial input dimension is small, the TCN expands the features into a higher-dimensional latent space before spatial diffusion, allowing the network to encode rich properties such as trends or local derivatives.

2.2.7 Optimization and structure learning

The model is trained to minimize the Mean Squared Error (MSE) between the predicted multi-horizon trajectories Y ̂ t + 1 : t + T out and the true observations Y t + 1 : t + T out. The objective function is defined in Equation 7 as:where ‖ ⋅ ‖ F denotes the Frobenius norm applied over the temporal horizon and the physical component dimensions.

L θ = 1 N ⋅ T out Y ̂ t + 1 : t + T out − Y t + 1 : t + T out F 2 , (7)

Since A is modeled as a continuous, differentiable parameter with a non-negativity constraint ( A i j ∈ R ≥ 0 ), we can learn the latent topology via standard gradient descent. Applying the chain rule, the gradient of the loss with respect to the structure A at layer l is given by Equation 8:where H temp represents the features processed by the TCN before spatial diffusion.

∂ L ∂ A = ∂ L ∂ H spat ⋅ ∂ H spat ∂ A = ∂ L ∂ H spat ⋅ H temp ⊤ , (8)

This mathematical derivation highlights a key insight: the gradient update for the edge weight A i j is proportional to the correlation between the error signal at node i and the temporal activity at node j. This mechanism effectively teaches the model to dynamically strengthen directed couplings that are useful for minimizing forecasting error, serving as an adaptive structural regularization.

3 Experimental setup

To validate the effectiveness and robustness of the proposed GNN-TCN architecture, we conducted a rigorous comparative analysis using historical geomagnetic data. This section details the data preparation, experimental design, and implementation parameters to facilitate architectural documentation and transparency. The implementation hyperparameters and training settings for the GNN-TCN architecture are summarized in Table 1.

TABLE 1. Summary of implementation hyperparameters and training settings for the proposed GNN–TCN architecture. The batch size and early-stopping patience were specifically tuned to stabilize convergence on the high-resolution MAGDAS dataset.

Specifically, to assess the model’s sensitivity to multi-scale temporal dependencies, we designed an experimental framework comprising four distinct model configurations. These configurations explore the trade-off between historical context and sampling granularity by varying the input window in physical time ( T i n ) and the downsampling factor ( D ), while maintaining a fixed forecasting horizon of 10 min. Consequently, the tensor sequence length fed into the network dynamically adapts to the resolution. The four scenarios evaluated are:

• Long-Context, Coarse-Grained (T i n = 600 min, D = 60s): Designed to capture macroscopic trends over a 10-h window with lower resolution (sequence length: 600 steps).
• Short-Context, Coarse-Grained (T i n = 60 min, D = 60s): Tests the model’s response to immediate history with low resolution (sequence length: 60 steps).
• Long-Context, Fine-Grained (T i n = 600 min, D = 30s): Provides a deep history with higher temporal density, challenging the model’s capacity to filter noise (sequence length: 1200 steps).
• Short-Context, Fine-Grained (T i n = 60 min, D = 30s): Focuses on high-frequency dynamics within the immediate past hour (sequence length: 120 steps).

3.1 Data preprocessing and spectral analysis

The dataset consists of high-resolution geomagnetic measurements from the MAGDAS station in Quito, Ecuador, covering the complete year of 2016. To strictly respect temporal causality, we avoided random shuffling and evaluated the models on a held-out chronological split.

3.1.1 Cleaning and chronological partitioning

Short gaps in the 1-min grid were linearly interpolated up to a maximum duration of 60 min, while remaining missing values were treated as break points. The cleaned series was partitioned into contiguous segments. We performed a chronological split over the ordered list of segments: Training Set (70%), Validation Set (15%, used for early stopping), and Test Set (Last 15%, strictly held out for final evaluation).

3.1.2 Normalization

Given the disparate dynamic ranges of the components (e.g., horizontal intensity H ≈ 27000 nT vs. declination variability), we applied Z-score normalization. The statistics ( μ train , σ train ) were computed exclusively using the training set and applied to transform the validation and test sets, preventing data leakage. All evaluation metrics are reported in physical units (nT) after applying the inverse transform.

3.2 Baseline models

We benchmark our proposed model against a comprehensive suite of classical and deep learning baselines to properly contextualize its performance. To ensure a fair comparison, all baselines were evaluated using the optimal temporal configuration established by our spectral analysis (T i n = 60 min, D = 60 s):

• Persistence (Naive Baseline): A standard zero-order hold benchmark assuming the field remains static over the prediction horizon (Y ̂ t + 1 : t + T out = Y t). This establishes the minimum physical threshold required to demonstrate that a model is learning dynamics rather than merely exploiting local inertia.
• Vector Autoregression (VAR): A classical, robust statistical baseline for multivariate time series forecasting that captures linear interdependencies among the time series components.
• LSTM (Long Short-Term Memory): A standard recurrent neural network baseline. To ensure a fair comparison and avoid initialization anomalies, the LSTM was systematically optimized via a hyperparameter grid search (evaluating hidden dimensions ∈ { 32,64,128 }, layers ∈ { 1,2 }, and learning rates ∈ { 1 e - 3,5 e - 4 }). The optimal configuration yielded a single-layer LSTM with 64 hidden units and a learning rate of 5 e - 4.

4 Results

We evaluated the predictive performance of the proposed GNN–TCN architecture against baseline models on a held-out chronological test split. Reported metrics include the Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and the Pearson Correlation Coefficient ( R ). Because the geomagnetic components ( H , D , Z ) operate on vastly different amplitude scales, the correlation was computed component-wise and then averaged to prevent artificial variance inflation. All errors are reported in physical units (nT) after inverse-transforming predictions using the training set statistics.

4.1 Ablation study: temporal context and resolution

To assess the model’s sensitivity to multi-scale temporal dependencies, we first evaluated the four distinct GNN–TCN configurations outlined in Section 3.

As shown in Table 2, the Short-Context, Coarse-Grained configuration yields the best overall performance among the deep learning setups. While the Naive Persistence baseline exhibits a marginally lower overall MAE (0.7858 nT vs. 0.8036 nT) due to short-term inertial artifacts at t = 1, the unoptimized standard LSTM baseline operates in a markedly higher error regime (MAE = 8.3013 nT). This initial discrepancy suggests that without proper structural regularization and hyperparameter tuning, standard sequential architectures struggle to disentangle the complex inter-component couplings from single-station data.

TABLE 2. Configuration ablation study on the test set (mean over H, D, Z) for a 10-min horizon. The unoptimized standard LSTM and Persistence baselines are compared against the four GNN–TCN spatial-temporal configurations.

4.2 Comprehensive baseline comparison and spatial ablation

Having established the Short-Context, Coarse-Grained setup as the optimal architecture, we expanded our evaluation against a broader suite of classical benchmarks and standard sequential deep learning baselines. To ensure absolute fairness, the LSTM baseline was systematically optimized via a hyperparameter grid search (single-layer LSTM, 64 hidden units, learning rate 5 e - 4), reducing its test error drastically to 1.9181 nT.

Crucially, to explicitly isolate the performance gain provided by the spatial graph regularization, we introduced a Pure TCN baseline. This model uses the exact same Dilated Temporal Convolution architecture as our proposal but entirely disables the Graph Learning and Graph Convolutional modules, processing each geomagnetic component independently.

As shown in Table 3, removing cross-component communication significantly degrades performance, increasing the MAE from 0.8036 nT (GNN–TCN) to 1.1259 nT (Pure TCN). This result primarily underscores the critical necessity of multivariate coupling in single-station modeling over isolated independent channel processing. While this ablation does not uniquely isolate the superiority of our learned graph topology over other dense channel-mixing multivariate structures, the high variance subsequently observed in the learned matrices (Section 4.5) suggests that treating this coupling through an adaptive graph formulation provides an effective, regularized feature-routing mechanism rather than a rigid, static mixing constraint. Overall, the proposed GNN–TCN successfully outperforms the Pure TCN baseline, the classical Vector Autoregression (VAR), and the fully optimized LSTM network.

TABLE 3. Comprehensive forecasting performance on the test set for a 10-min horizon. The optimal GNN–TCN is compared against the Pure TCN ablation, as well as optimized classical and deep learning baselines.

It is worth noting an apparent metric discrepancy in Table 3: the Persistence baseline achieves a marginally lower MAE (0.7860 nT) than the proposed GNN–TCN (0.8036 nT), yet exhibits a substantially lower correlation coefficient (R = 0.8143 vs. R = 0.9950). This discrepancy arises from the fundamentally different aspects of forecast quality that each metric captures. The Persistence model is a trivial baseline that propagates the last observed value forward across the entire prediction horizon. This strategy incurs low absolute errors during the predominant periods of geomagnetic quiescence, where the signal changes slowly and the last observation remains a locally valid approximation. However, this same strategy fails to track the temporal dynamics of the signal, producing a systematic phase lag that severely penalizes its correlation with the true trajectory. The correlation coefficient, being sensitive to the shape and phase alignment of the predicted signal rather than its absolute magnitude, therefore reveals that Persistence offers a structurally poor representation of the geomagnetic field evolution. In contrast, the GNN–TCN actively models the joint non-linear evolution of the ( H , D , Z ) components, achieving near-perfect correlation while incurring only a marginal increase in mean absolute error. These two metrics are therefore complementary: MAE quantifies the average magnitude of deviations, while R quantifies the fidelity of the predicted temporal structure. Taken together, they confirm that the GNN–TCN provides a structurally superior forecast, particularly during the dynamic transients that are most critical for operational space weather applications.

To ensure statistical significance and rule out initialization anomalies, we conducted a robustness analysis by retraining the optimal GNN–TCN architecture across 10 distinct random seeds. The model achieved a robust average MAE = 0.8178 ± 0.1383 nT, RMSE = 1.1866 ± 0.1281 nT, and R = 0.8084 ± 0.0117 across all runs (with our representative seed achieving the reported 0.8036 nT). This scale difference in the correlation coefficient ( R ) arises from a methodological distinction in metric aggregation: the representative value of 0.9950 in Table 3 is a global Pearson correlation, computed by flattening the complete continuous year-long testing trajectory, where large macro-scale diurnal variations dominate the variance and guide the linear trend. Conversely, the multi-seed robustness script calculates the correlation locally window-by-window across individual multi-horizon segments and averages them. This windowed approach serves as a far more stringent metric for micro-scale perturbations but naturally shifts the mean to a lower baseline scale. Both scales, however, independently confirm stable and highly consistent phase tracking across initializations, as evidenced by the tight multi-seed standard deviation.

4.3 Error drift and temporal stability

To rigorously assess the predictive stability of the models, we analyzed the error degradation rate (Error Drift) across the 10-min forecasting horizon. As depicted in Figure 4 and supported by numerical analysis, classical and deep learning architectures exhibit fundamentally different behaviors when projecting future geomagnetic states.

The standard Persistence model achieves the lowest absolute error at t = 1 minute (MAE = 0.2016 nT). This is an expected artifact in high-resolution space weather forecasting, as the geomagnetic field exhibits extreme localized autocorrelation; merely carrying forward the last known observation serves as a mathematically strong, albeit physically uninformative, short-term baseline. However, as the horizon expands, the limitations of non-dynamic models become apparent. Persistence and VAR exhibit rapid degradation rates, accumulating + 1.0936 nT and + 1.0551 nT of error by t = 10, respectively.

Conversely, the proposed GNN–TCN demonstrates superior temporal stability, registering the lowest overall error degradation rate (+ 0.8451 nT). Rather than relying on immediate local inertia, the GNN–TCN learns the underlying non-linear coupling dynamics, allowing its error curve to grow at a significantly slower and more controlled rate. By integrating Dilated Temporal Convolutions (TCN) alongside our learned graph structure, our architecture maintains the necessary localized receptive fields to capture sudden transients without accumulating linear drift, confirming its suitability for continuous operational forecasting.

4.4 Model robustness under disturbed geomagnetic conditions

To evaluate the operational readiness of the proposed model during severe space weather events, a targeted stress test was conducted by segregating the test set into two distinct subsets based on the temporal variance of the target window: Quiet Conditions (bottom 95% of the data) and Storm Conditions (top 5% of maximum variance, corresponding to sudden storm commencements and active phases). Figure 5 details the comparative performance.

Under highly predictable Quiet Conditions, classical benchmarks leverage local signal inertia, with the zero-order hold Persistence baseline yielding a low error (MAE = 0.6931 nT), slightly outperforming the statistical VAR (0.7794 nT) and the proposed GNN–TCN (0.8147 nT). This behavior is physically consistent with an unperturbed magnetosphere, where vector fluctuations remain minimal and stationary.

However, during Storm Conditions, the system shifts into highly non-linear regimes where inertial assumptions fail. The baseline models exhibit severe degradation; the Persistence error amplifies by more than 3.6-fold, jumping to 2.5496 nT, while the thoroughly optimized LSTM approach collapses entirely, accumulating an extensive error of 3.4546 nT due to its inherent phase lag during high-amplitude transients.

In sharp contrast, the proposed GNN–TCN architecture demonstrates notable resilience under disturbed conditions, achieving the lowest absolute error among all evaluated configurations (MAE = 2.2670 nT). By regularizing space weather predictions through a learned dynamic graph structure and processing multi-scale history with dilated inception kernels, the GNN–TCN filters high-frequency stochastic noise while tracking the macro-scale structural transitions of the storm. This consistent mitigation of error spikes under intense turbulence suggests that our topological approach represents a promising candidate for integration into real-time space weather alerting pipelines.

4.5 Topological stability and adaptive regularization

The learned adjacency matrices, A, provide a transparent view into the spatial dependencies inferred by the GNN–TCN architecture. Figure 6 presents the mean ( μ A ) and standard deviation ( σ A ) of the adjacency matrices across 10 independent training runs with different random initializations.

The mean adjacency matrix (Figure 6, left panel) reveals a learned topological prior, characterized by a strong coupling between the H and D components, consistent with their physical interplay in the horizontal plane. Conversely, the standard deviation heatmap (Figure 6, right panel) exhibits high variance (σ A ≈ 0.4076 for the highest connection), resulting in a mean Coefficient of Variation of 102.43% across the runs.

This high variance is a pivotal finding: it explicitly proves that the GNN–TCN does not collapse into a rigid, uniquely identifiable static physical law. Instead, the graph structure functions as an adaptive structural regularizer, allowing the model to dynamically route temporal features across the ( H , D , Z ) nodes. This stochastic adaptability provides a topological inductive bias that stabilizes the forecasting process and minimizes generalization error, allowing our architecture to navigate non-linear coupling more effectively than fixed-structure models.

5 Discussion

The experimental results provide crucial insights into why classical statistical models and standard deep sequence architectures struggle in the single-station regime, and how explicit topological structure improves multi-horizon geomagnetic predictions. The clear separation in error regimes across the forecasting horizon highlights the importance of tailored structural inductive biases. Unstructured architectures, such as the optimized LSTM baseline, treat the multivariate input [ H , D , Z ] as a flat feature space, forcing them to learn spatiotemporal correlations implicitly without directional constraints. In a highly chaotic system like the local geomagnetic field, this lack of structured coupling can lead to a rapid accumulation of systematic errors and vulnerability to phase lags during high-amplitude transients. While a zero-order hold baseline like Naive Persistence yields lower absolute distance errors (MAE) during dominant quiet-time intervals by exploiting extreme localized field autocorrelation, it lacks the dynamic capacity to track the field’s underlying physical evolution, resulting in severe phase degradation over expanding horizons. In contrast, the GNN–TCN architecture explicitly separates temporal feature extraction from inter-component dependency modeling, ensuring that the dimensions exchange information through a constrained routing mechanism that mitigates multi-step error drift without assuming rigid physical identities.

Crucially, our analysis of the learned adjacency matrices across multiple random initializations revealed high variance, demonstrating that the model does not converge to a single, unique physical interaction graph. This finding supports the interpretation that the performance gains of the GNN–TCN stem from treating multivariate coupling not through a rigid, dense channel-mixing layer, but through the graph structure acting as an adaptive structural regularizer. Rather than discovering a fixed, immutable physical law, this flexible topological routing induces a data-driven regularization that stabilizes chaotic single-station dynamics.

From a geophysical perspective, this adaptive regularization is highly consistent with the nature of single-station observations. Without the explicit inclusion of global upstream drivers—such as solar wind parameters or the Interplanetary Magnetic Field (IMF)—a model cannot fully map the global dynamics of the solar-terrestrial coupling. Instead, the inferred couplings loosely reflect local transient phenomena. For instance, the strong interaction generally observed between the H and D components correlates with the horizontal plane dynamics driven by the Equatorial Electrojet (EEJ), which produces a pronounced daytime enhancement of the magnetic field near the magnetic equator. The vertical component Z, meanwhile, adjusts dynamically based on local induced currents within the Earth’s conducting interior. Within this context, the learned graph topology provides an effective, data-driven framework to operationalize these multivariate couplings transiently, allowing the network to filter stochastic noise while locking onto deterministic low-frequency trends.

Beyond short-term prediction accuracy, this framework holds substantial operational relevance for space weather monitoring. In practical applications, the most critical quantity is the temporal rate of change of the magnetic field ( d B / d t ), which drives the magnitude of geoelectric fields induced in the ground. Crucially, while a static Naive Persistence baseline can minimize absolute distance metrics during quiet periods, its zero-order hold assumption inherently yields a zero rate of change ( Δ B / Δ t = 0 ) between consecutive predicted steps, rendering it operationally blind to the rapid transitions that trigger Geomagnetically Induced Currents (GIC). These sudden variations pose severe systemic risks to power transmission networks and pipelines. By effectively modeling the joint non-linear evolution of ( H , D , Z ), and demonstrating robust predictive resilience during the top 5% of highly disturbed storm transients, the proposed GNN–TCN architecture provides an operationally viable mathematical building block for real-time risk assessment and alerting systems.

5.1 Limitations and future work

As a proof-of-concept study, we acknowledge that the reliance on single-station data without external drivers inherently constrains the model’s interpretation of global magnetospheric drivers. Furthermore, while the learned graph topology over ( H , D , Z ) effectively operationalizes these cross-component couplings, it remains an architectural design choice—acting as an adaptive structural regularizer—rather than a strongly grounded physical derivation. We also note as an architectural limitation that this baseline comparison was restricted to standard sequential and statistical benchmarks (Persistence, VAR, LSTM), leaving the systematic evaluation against large-scale Spatio-Temporal Graph Neural Networks (STGNNs) or heavy multi-station temporal transformers for future framework extensions. Additionally, reproducibility remains constrained by operational factors, as the high-resolution MAGDAS observatory data is subject to institutional distribution restrictions, and open-source code repositories are not publicly available at this proof-of-concept stage. Consequently, this framework is presented strictly as a localized pilot framework. Extending this approach to multi-station global networks, integrating upstream solar wind parameters, and validating across multiple solar cycles will be necessary to transition this proof-of-concept into a broadly generalizable forecasting infrastructure.

6 Conclusion

This paper introduces a topologically informed proof-of-concept approach for short-term, single-station geomagnetic forecasting. By integrating a learnable directed graph over the ( H , D , Z ) vector components with a Temporal Convolutional Network (TCN) based on causal dilated convolutions, the proposed GNN–TCN explicitly captures structured cross-component interactions while maintaining highly stable multi-horizon predictions.

The primary findings and contributions of this study are summarized as follows:

• Adaptive Structural Regularization: We demonstrated that representing the ( H , D , Z ) components as nodes in a learned directed graph enables the architecture to capture asymmetric inter-component interactions. Extensive multi-seed evaluations proved that this adjacency matrix acts as a highly effective, adaptable structural regularizer rather than discovering a rigid, static physical law, ultimately mitigating multi-step error drift.
• Fidelity Over Baseline Benchmarks: On a held-out chronological test split, the GNN–TCN architecture established a robust predictive framework, achieving a global Mean Absolute Error (MAE) of 0.8036 nT and exceptional phase alignment ( R = 0.9950 ). While the Naive Persistence baseline capitalizes on localized field autocorrelation during dominant quiet-time intervals to minimize absolute distance metrics, our topological model successfully handles structural transitions where zero-order assumptions fail, offering competitive tracking capabilities against classical Vector Autoregression (VAR) and thoroughly optimized LSTM baselines.
• Resilience in Disturbed Conditions: Unlike standard sequential baselines that collapsed under the non-linear volatility of extreme space weather, the GNN–TCN maintained predictive stability during a targeted stress test, gracefully controlling error degradation and achieving the lowest absolute forecasting error (MAE = 2.2670 nT) during the top 5% most severe local geomagnetic storm events.

While this framework establishes a robust baseline for single-station geomagnetic modeling, its current scope intentionally isolates local topology from global phenomena. By relying strictly on local magnetic perturbations without external solar-terrestrial drivers (e.g., solar wind parameters or IMF), the model infers proxy dynamics rather than the full physical solar-magnetosphere coupling. Furthermore, empirical validation was restricted to a single equatorial station over a one-year period. Future work will naturally extend this architectural framework to multi-station networks, enabling the capture of spatially distributed geomagnetic dynamics. Integrating external upstream drivers and evaluating cross-year robustness spanning different phases of the solar cycle will further bridge the gap between data-driven topological learning and operational space weather forecasting.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Yuichi Otsuka, Nagoya University, Japan

Reviewed by: Mohamed Freeshah, Wuhan University, China

Aiganym Sarsembayeva, Al-Farabi Kazakh National University, Kazakhstan

Henok Tenaw Moges, University of Cape Town, South Africa