Domain-Informed Graph Neural Networks for Climate Factor Forecasting to Support Sustainable Crop Management

1 Introduction

Climate factors, such as temperature, precipitation and humidity, are climatic drivers that directly influence crop phenology, soil moisture dynamics, and pest & disease pressures. Therefore, predicting changes in these factors is irreplaceable for efficient crop management and food security planning [Nouri and Veysi, 2025]. However, achieving high-resolution and robust predictions remains challenging because local climatic dynamics are inherently influenced by surrounding regions. Without explicitly considering the spatial dependencies among neighboring areas, predictions often fail to capture propagating weather patterns and regional interactions, especially at large spatial scales. Moreover, climate regimes vary substantially across regions, exhibiting non-uniform feature distributions and heterogeneous temporal patterns, which makes it difficult for a single model to generalize effectively across diverse agro-climatic zones.

Traditional climate forecasting has long relied on numerical weather prediction (NWP) and climate models grounded in atmospheric dynamics [Kimura, 2002; Leutbecher and Palmer, 2008; Zhao et al., 2025; Giorgi et al., 2023]. These methods have achieved strong performance in large-scale weather and climate simulation by explicitly modeling physical processes. However, their effectiveness often depends on complex physical parameterizations, high-quality boundary conditions, and substantial computational resources, which may limit their flexibility for regional, station-level, and application-oriented agro-climatic forecasting.

To complement physics-based approaches, data-driven models have been increasingly explored for meteorological and climate prediction. Classical machine learning methods such as SVMs and random forests have been used to model nonlinear relationships among meteorological variables [Das and Ghosh, 2018; He et al., 2021b]. Deep learning models, including CNNs, RNNs, ConvLSTMs, and recent spatiotemporal architectures, further improve the ability to capture temporal evolution and spatial patterns from historical observations [Haggag et al., 2021; Chen et al., 2022; Gong et al., 2024; Utku and Akyol, 2025]. Despite these advances, many existing methods still rely mainly on statistical correlations and may insufficiently encode domain-specific agro-climatic structure.

Recently, graph neural networks (GNNs) have shown strong potential in climate-related forecasting by representing stations, grids, or regions as nodes and modeling spatial dependencies through graph propagation [Humphries et al., 2026; Yousaf et al., 2025; Zhang et al., 2026; Li et al., 2025]. Such models are particularly suitable for multi-station forecasting because they can capture non-Euclidean spatial interactions beyond regular grid assumptions. Nevertheless, most existing graph-based methods construct spatial graphs from geographical distance, statistical similarity, or learnable adjacency matrices, while the integration of agro-climatic domain knowledge into graph construction remains underexplored.

Compared with NWP and traditional machine learning methods, deep learning methods-especially graph-based models-exhibit improved performance in capturing nonlinear spatiotemporal correlations of univariate data. Nevertheless, the prediction of climate data via deep learning still encounters the following challenges:

1. Existing GNN-based methods construct graphs using geographic proximity, statistical similarity, or fully learnable adjacency. Such topologies are often misaligned with agro-climatic regimes.
2. Existing models simply encode all components in a single representation, which may cause representation redundancy.
3. Many methods capture only one type of dependency (static or dynamic) or yield opaque learned edges without clear agro-climatic meaning.

Based on the above challenges, we propose DoIGNN, a Domain-Informed Graph Neural Network, a framework that explicitly integrates agro-climatic domain knowledge into graph-based spatiotemporal modeling. In this approach, climate observation stations are systematically grouped into Agro-Climatic Homogeneous Zones (ACHZs) using long-term meteorological statistics and geospatial attributes. Specifically, intra-zone connections capture localized climatic coherence, while inter-zone links selectively model broader regional interactions and teleconnections. Furthermore, the spatiotemporal decomposition module is introduced to factorizes the multivariate climate tensor into a low-rank set of global temporal bases and station-specific loadings, and then derives a station-level global component that summarizes common spatiotemporal trends for target forecasting. DoIGNN jointly performs prediction over the ACHZ-based and learned graphs, thereby uncovering both cross-variable and cross-regional dependencies. The environmental-resource science domain experts from Ningbo Natural Resources and Planning Big Data Center codeveloped the ACHZ construction by specifying agronomically meaningful indicators and thresholding rules. In addition, they provided domain-relevant criteria for more comprehensive and accurate model evaluation.

The main contributions of this article are as follows.

1. DoIGNN partitions stations into ACHZs and constructs a ACHZ-guided hierarchical graph that encourages within-zone edges while selectively allowing cross-zone connections.
2. By spatiotemporal decomposition, DoIGNN exploits shared climate dynamics while avoiding excessive reliance on raw auxiliary variables. The temporal regularization further preserves temporal structure, leading to more stable learning and improved robustness.
3. DoIGNN performs forecasting on real-word climate datasets validate the effectiveness of our model. Compared with several baseline methods under different prediction horizons, our method demonstrates superior performance.

2 Preliminaries and Problem Statement

2.1 Agro-Climatic Homogeneous Zones (ACHZs)

Agro-Climatic Homogeneous Zones (ACHZs) partition monitoring stations into groups sharing similar long-term agroclimatic regimes (e.g., thermal resources, moisture availability, and rainfall seasonality). Unlike purely data-driven clustering, ACHZs are defined by rule-based agro-climatic indicators computed from multi-decadal daily meteorological records, improving interpretability and aligning the graph structure with agro-ecological principles [Fischer et al., 2012].

Given the annual accumulated Growing Degree Days GDDi, the Aridity Index AIi, and the Seasonality Index SIi, each station i can be assigned to an ACHZ by discretizing these indicators into ordered classes with agronomically meaningful thresholds:

where ϕth, ϕmo, ϕse are rule-based binning functions (e.g., low/medium/high heat, arid-humid regimes, weak-strong seasonality), and ⊕ denotes concatenation into a unique ACHZ identifier. This yields zones that are directly interpretable in terms of agricultural climate constraints, enabling principled spatial regularization in subsequent graph modeling.

2.2 Problem Statement

Consider a multivariate spatiotemporal climate dataset X ∈ RT×N×C , where T denotes the number of time steps (days), N is the number of stations, and C is the number of climate factors. Given a target channel, the goal is to forecast its future values Y ∈ Rf×N over a future horizon of f = 12 days using the historical h = 12 days of all variables.

Formally, for each station n ∈ {1, . . . ,N} and each reference time t satisfying L ≤ t ≤ T -H, we define the input history tensor and prediction target as:

We aim to learn a forecasting function fθ parameterized by θ such that:

where G denotes a set of graphs describing inter-station dependencies.

3 Methods

The DoIGNN architecture mainly consists of ACHZs Construction, graph generation and prediction head module. The schematic diagram of DoIGNN can be seen in Figure 1.

3.1 ACHZs Construction

We construct an ACHZ-guided graph to inject domain knowledge into spatial dependency modeling. Each station is first assigned a zone label using rule-based discretization of longterm agro-climatic indicators, and then a weighted adjacency matrix is derived by combining within-zone consistency with cross-zone similarity.

Let D(y) i denote the set of days in year y at station i. Using daily mean temperature Ti(d) (◦C), daily minimum/maximum temperatures T min i (d), T max i (d) (◦C), daily precipitation Pi(d) (mm), and daily incoming solar radiation Rs,i(d) (MJ m-2 day-1 ), we compute the following longterm descriptors.

Station-Wise Agro-Climatic Descriptors

For station i, we compute a climatological descriptor vector

where Ti is the multi-year mean daily temperature, GDDi(Tb) is the multi-year mean annual GDD with base temperature Tb, FFPi is the mean frost-free period, HDi is the mean number of heat-stress days (e.g., T max ≥ 35◦C).

Thermal Resources via GDD

For a crop-specific base temperature Tb, the annual accumulated GDD is

We use the multi-year climatological mean GDDi = 1 Y ∑y GDD(y) i as a proxy of heat availability, which is widely used to quantify crop development and phenological progress.

Moisture Availability via AI

Following a commonly adopted definition in dryland and agro-environmental studies, we define the aridity index as

where Pi is the multi-year mean annual precipitation (mm year-1 ) and PETi is the multi-year mean annual potential evapotranspiration (mm year-1 ). The AI = P/PET formulation is standard in global aridity characterization [Middleton et al., 1992].

When full meteorological inputs for FAO-56 PenmanMonteith are unavailable, we estimate PET using the Makkink reference evapotranspiration estimator based on solar radiation and pressure-related psychrometric variables, and then aggregate daily ET0 to annual PET [Allen et al., 1998].

Rainfall Seasonality via SI

To characterize intra-annual rainfall concentration, we compute a seasonality index based on monthly precipitation climatology {P i,m}12 m=1 :

where P i and PETi are multi-year mean annual precipitation and potential evapotranspiration, and {Pi,m}12 m=1 is monthly precipitation climatology. SI measures deviation from uniformly distributed rainfall across months and is commonly used to describe rainfall seasonality patterns [Walsh and Lawler, 1981].

3.2 Spatiotemporal Decomposition

Given a 3-D tensor Y ∈ RT×N×C , the decomposition module separates Y into a shared global component and a local residual:

We model Y(g) via a rank-R CANDECOMP/PARAFAC (CP) factorization [Harshman, 1970; Carroll and Chang, 1970; Kolda and Bader, 2009]:

where A ∈ RN×R captures spatial loadings, B ∈ RC×R captures variable loadings, and X ∈ RR×T denotes global temporal factors shared across all stations and variables. Equivalently, Y(g) = ∑Rr=1 ar ◦ br ◦ xr with ◦ being the outer product. The local term is defined as the residual Y(l) = Y - Y(g) .

To encourage temporally consistent dynamics in X, we introduce a one-step predictive regularizer using a Temporal Convolutional Network (TCN) ϕ(·) [Bai et al., 2018]. For a window length w, we penalize deviations from the one-step prediction:

Finally, the factors are estimated by minimizing a regularized reconstruction objective:

where λ, µ > 0 control weight decay and temporal smoothness, respectively. In practice, we solve ((11)) by alternating updates: fixing (B,X) to update A, fixing (A,X) to update B, and fixing (A,B) to update X, while periodically applying the TCN-based temporal step.

For downstream forecasting, we derive a station-wise global representation as auxiliary information by projecting temporal factors through spatial loadings:

To avoid information leakage, the decomposition is fitted using only the training period.

3.3 Pre-trained Encoder

Our encoder follows the masked autoencoder (MAE) paradigm [He et al., 2021a]. Given the input sequences, we perform patch-wise masking and representation learning for each station. For node i, let Xi ∈ RT denote the input sequence, we split it into P non-overlapping patches of length L, where the j-th patch is Xi,j ∈ RL. Then, a high ratio of patches is randomly masked and fed into the corresponding encoder. The encoder is composed of four Transformer blocks to capture temporal patterns, producing patch embeddings Hi j ∈ Rd for local inputs.

Before decoding, masked positions are filled with trainable mask tokens. The decoder contains one Transformer block followed by an MLP head, which reconstructs the j-th patch as ˆXij ∈ RL from Hi j . We optimize the reconstruction error only on masked patches using RMSE:

where M is the number of masked patches, Mm is the corresponding mask matrix, ∥ · ∥F denotes the Frobenius norm, and sum(·) sums all entries of a matrix. After training, the trained encoder can be used to model the overall understanding of the time series, as the autoencoder masks the majority of the input sequence to encourage the learning of high-level semantics.

3.4 Graph Generation

To fully model spatiotemporal dependencies in climate data, DoIGNN constructs (i) a dynamic-static adjacency from node embeddings and raw sequences, and (ii) an ACHZ-guided weighted adjacency. The flowchart is shown in Figure 2.

Dynamic-Static Adjacency

For node i, we extract a static representation from the input sequence Xi via a lightweight CNN-MLP encoder:

where each Conv(·) denotes a convolutional block with convolution, ReLU, and batch normalization, and FC(·) denotes a fully-connected block with ReLU and batch normalization. We then fuse the static feature Gi with the masked embedding Hi to obtain a dynamic-static representation:

We aim to learn a sparse adjacency by modeling edge existence between node pairs. Given Z i and Z j , we compute the Bernoulli logits and probabilities as:

where Θij ∈ R2 denotes unnormalized logits and Θ′ij is the corresponding normalized Bernoulli distribution. To enable end-to-end optimization with discrete sampling, we employ the Gumbel-Softmax reparameterization [Jang et al., 2022]:

where g ∈ R2 is sampled i.i.d. from Gumbel(0, 1) and τ is the temperature. To regularize the learned structure and control sparsity, we introduce a kNN prior graph Ak and penalize deviations using a cross-entropy term:

where Ak specifies the kNN connectivity.

ACHZ-Guided Weighted Adjacency

Let C = [c1, . . . , cN]⊤ be the descriptor matrix. We first standardize each feature dimension:

and define an agro-climatic similarity kernel

where σ controls the kernel bandwidth. The ACHZ-guided weighted adjacency AACHZ is

where λ ∈ (0, 1) down-weights cross-zone connections to emphasize within-zone homogeneity while still allowing weak long-range interactions.

3.5 Prediction Head

We adopt a downstream GNN, an ACHZ-guided Graph WaveNet, to performmulti-step forecasting. To jointly leverage the learned dynamic-static ACHZ-guided weighted adjacency, we treat Ads and AACHZ as undirected graph and obtain forward & backward transitions.

ACHZ-guided Graph WaveNet further employs dilated causal convolutions to capture temporal dependencies. In the fusion layer, we combine the latent representation with the last patch embedding HP to obtain the final prediction:

where Xtarget ∈ RN×T is the target factor. || denotes concatenation, MLP consists of two fully-connected layers with ReLU, and Conv is a 1×1 convolution with ReLU. The output ˆY represents the f-step-ahead forecasts for n nodes.

Finally, we optimize the prediction model using MAE together with the graph regularizer in:

where Y is the ground truth, and λgraph controls the strength of structure regularization. An overview of the full training pipeline is provided in Algorithm 1.

4 Experiments

4.1 Dataset

We use the China Meteorological Assimilation Driving Datasets for the SWAT model-Long time series (CMADS-L V1.0), obtained from the official CMADS portal. 1 CMADS

1 https://www.cmads.org/

Algorithm 1: DoIGNN

Input: Past observations: X ∈ RT×N×C . Output: Future forecasts: ˆY ∈ Rf×N . /* Spatiotemporal Decomposition */ 1 for i1 alternating iterations do 2 Fix TCN ϕ(·),A,B update X using Equation (11). 3 Fix TCN ϕ(·),X and update A,B using Equation (11). 4 if mod iteration then 5 Fix A, B, X and update TCN ϕ(·) using Equation (10). 6 end 7 end 8 Obtain station-wise global representation Xaux by Equation (12). /* Pre-trained Encoder */ 9 Split Xinto P patches and mask them with probability p. 10 for i2 epochs do 11 Encode X to obtain H, complement masked positions with trainable patches, and decode to yield ˆX. 12 Update encoder & decoder parameters using Equation (13). 13 end 14 Obtain pretrained node embeddings H. /* Prediction Head */ 15 Compute ACHZ through Equation (4) - (7) 16 Construct ACHZ-guided weighted adjacency AACHZ by Equation (19) - (21). 17 for i3 epochs do 18 Generate dynamic-static adjacency Ads by Equation (14)-(18). /* Prediction module */ 19 Compute transition matrices and perform dilated causal convolutions by Equation (22) to obtain Hgw. 20 Fuse HP and Hgw via (22) to produce ˆY. 21 Update prediction head parameters using (23). 22 end

is a meteorological forcing dataset designed for hydrological and agro-environmental applications, produced by assimilating multi-source observations and numerical background fields [Meng et al., 2018]. CMADS V1.0 provides daily meteorological fields over East Asia (approximately 0◦ - 65◦N, 60◦ -160◦E) at 1/3◦ spatial resolution, organized as a 195×300 grid. Climate factors includes: surface atmospheric pressure p (hPa), daily mean temperature ¯ T (◦C), daily maximum/minimum temperatures T max , T min (◦C), 24-hour accumulated precipitation P (mm), relative humidity RH (%), solar radiation Rs (MJm-2 ), specific humidity q (gkg-1 ), and wind speed u (m s-1 ).

In this study, the forecasting target is the daily mean temperature, solar radiation and specific humidity. We subsample N = 130 stations on a regular grid by selecting one site every 15 latitude indices and one site every 30 longitude indices, yielding a 13 × 10 station set.

4.2 Evaluation Criteria and Experiment Settings

Beyond point-wise error metrics, we adopt two metrics that are widely used to assess model performance in environmental science and meteorology: Nash-Sutcliffe Efficiency (NSE) [Nash and Sutcliffe, 1970] and the Kling-Gupta Efficiency (KGE) [Gupta et al., 2009].

Let yt(i) and ˆ yt(i) denote the observed and predicted values at station i and time index t over an evaluation window T , and let ¯ y(i) = 1 |T | ∑t∈T yt(i) be the sample mean. The NSE is defined as

NSE quantifies the relative improvement over a mean-value baseline [Nash and Sutcliffe, 1970].

To provide a more diagnostic assessment that disentangles error sources, we further report KGE, which jointly evaluates linear association, bias, and variability [Gupta et al., 2009]:

where r(i) is the Pearson correlation between {yt(i)} and {ˆ yt(i)}, β(i) = µˆ y(i)/µy(i) measures the mean bias ratio, and γ(i) = σˆ y(i)/σy(i) measures the variability ratio, with

µ(·) and σ(·) denoting the sample mean and standard deviation over T . Similar to NSE, values closer to 1 imply better overall agreement [Gupta et al., 2009].

In spatiotemporal decomposition, the learning rate is 0.0005. The temporal regularization parameter λGLD is 0.1. TCN ϕ(·),A,B, X are optimized in 50 iterations and alternate training fashion for every 5 iterations. In embedding module, mask probability p is set as 75%. Embedding and prediction module are both trained for i2 = i3 = 100 times. The learning rate for embedding is 0.0005 without decay, and for prediction is 0.005 with decay. The softmax temperature parameter τ is 0.5. Finally, The graph regularization parameter λgraph is 0.1. Details for ACHZs can be seen in Table 3. Relevant code is available at the project repository.

4.3 Forecasting Results

Figure 3 visualizes the constructed ACHZs, which provide an interpretable partition of stations with similar long-term agroclimatic regimes. Table 1 reports the macro-averaged multihorizon forecasting results over all stations, where the best scores are highlighted in bold. Overall, DoIGNN achieved the most consistent performance across horizons and variables, especially on skill-oriented metrics (NSE and KGE), indicating improved predictive reliability beyond point-wise errors.

Specifically, in average temperature dataset, DoIGNN attains the best overall skill for temperature forecasting, yielding the highest NSE and KGE at all horizons. Notably, some baselines like WaveNet can obtain relatively small MAE but exhibit substantially degraded NSE and KGE. This discrepancy suggests that those models without graph information may fit the central tendency while failing to reproduce the temporal variability and correlation structure. GWNet, on the other hand, markedly improved prediction errors by integrating GCN rather than TCN for spatial relationship analysis. This underscored the significance of graph relationships from domain knowledge in augmenting predictive performance of climate factor.

Finally, Table 2 demonstrates the ablation study results. Removing the ACHZ-guided adjacency (AACHZ ) causes a clear degradation in NSE/KGE, indicating that domainconsistent spatial constraints are critical for preserving forecast skill and preventing spurious cross-regime message passing. Removing the spatiotemporal decomposition (STD) further deteriorates performance, confirming that explicitly extracting global shared dynamics improves stability and generalization. When both STD and AACHZ are removed, the model collapses to the WaveNet-style backbone and exhibits the same behavior as the corresponding baseline, demonstrating that the performance gains of DoIGNN primarily come from the proposed domain-informed graph prior and the global-local decomposition rather than architectural scaling alone.

5 Discussion

The proposed forecasting model can be deployed as a decision-support component for crop management. Given multi-station predictions of climate factors over the forthcoming horizon, the system can derive risk-aware indicators at the zone and field scales such as heat & cold stress likelihood, rainfall deficit & excess, and humidity-driven disease favorability. These indicators enable practitioners to (i) schedule irrigation and drainage by anticipating short-term water stress, (ii) adjust sowing, fertilization, and harvest windows to avoid adverse temperature/rainfall regimes, and (iii) trigger early warnings for pest and fungal disease management under highrisk humidity-temperature patterns. By aggregating forecasts within agro-climatic homogeneous zones, the framework also supports coherent regional planning.

6 Conclusion

This paper presented DoIGNN, a domain-informed graph neural network for climate factor forecasting. DoIGNN injects agro-climatic knowledge via an ACHZ-guided adjacency, and further employs a temporally-regularized globallocal spatiotemporal decomposition to obtain compact shared dynamics as auxiliary information. Built upon the ACHZ constraint and a learned dynamic-static dependency, the proposed prediction model jointly captures cross-region and cross-variable interactions. Experiments on CMADS demonstrate consistent improvements over strong baselines across horizons. In the future, we will focus on lightweight and scalable modeling for long-horizon forecasting over large-scale climate archives to reduce computational and storage costs without sacrificing predictive skill.

Ethical Statement

There are no ethical issues.

Acknowledgements

This work is supported by National Key R&D Program of China (Grant No. 2021YFC2101100), National Natural Science Foundation of China (62073288, 12075212), and Zhejiang University NGICS Platform, and their supports are thereby acknowledged.