M$^2$Weather: A Benchmark for Joint Multi-Station and Multi-Variable Weather Forecasting

1 INTRODUCTION

• Existing studies often separate station and variable interactions and use differing datasets and settings, leaving their individual and joint contributions without systematic comparison.The benchmark is motivated by the need for high-quality observations and physically coupled weather variables to assess these relationships reliably.Source: §1 · p. 2
• M2Weather provides 2,809 high-quality stations, five physically coupled variables, three spatial scales, and unified forecasting protocols.The France, Europe, and Global settings support fair comparisons across network extents.Source: §1 · p. 2
• The lightweight, plug-and-play adapter adds missing station or variable interactions to a trained model without retraining, enabling controlled analysis of each dimension.The adapter complements relationships missing from a baseline.Source: §1 · p. 2
• Across the three datasets, MSMV methods achieve the best performance, and every adapted backbone improves consistently with interaction completion.The benchmark evaluates 16 baselines; completing all station–variable interactions lets some non-MSMV methods achieve state-of-the-art performance.Source: §1 · pp. 2–3

2 TAXONOMY OF STATION WEATHER FORECASTING

• The taxonomy classifies methods by explicit cross-station and cross-variable interactions; shared parameters or pre-training alone do not imply either interaction.Figure 2 also outlines the benchmark workflow: station stratification, wrappers preserving backbone interaction structures, and adapter analysis on frozen backbones.Source: §2 · p. 3, Figure 2 · p. 4

2.1 SINGLE-STATION SINGLE-VARIABLE FORECASTING

• SSSV methods model individual series without explicitly coupling station or variable observations, providing temporal-modeling baselines.Examples span channel-independent models and foundation models that treat data as single-series sequences.Source: §2.1 · p. 3

2.2 SINGLE-STATION MULTI-VARIABLE FORECASTING

• SSMV methods exploit dependencies among quantities within one station but do not explicitly represent station-level spatial structure.Examples include attention across variables, channel relationships, and any-variate or multivariate forecasting models.Source: §2.2 · p. 3

2.3 MULTI-STATION SINGLE-VARIABLE FORECASTING

• MSSV methods exploit relationships among stations for one weather variable, spanning explicit spatial interaction and lightweight spatial or temporal knowledge.Examples include multi-scale spatial correlation, subgraph attention, field reconstruction, and prompt- or coordinate-based approaches.Source: §2.3 · p. 3

2.4 MULTI-STATION MULTI-VARIABLE FORECASTING

• MSMV methods jointly model within-station variable coupling and cross-station dependencies, a comparatively underexplored approach.Examples include adaptive variable and station graphs, cross-regional dynamics with higher-order variable interactions, and station-derived surface fields.Source: §2.4 · p. 3

3.1 PROBLEM DEFINITION

• The forecasting task maps historical observations across N stations and V variables to future targets using a model fθ.The input and target tensors have shapes L×N×V and H×N×V, respectively.Source: §3.1 · p. 4

3.2.1 DATA COLLECTION

• The source product supplies land-surface station observations; M2Weather extracts five physically coupled variables from 2017–2021, yielding 17,173 candidate stations and 43,824 hourly timestamps.The extracted quantities are air temperature, dew-point temperature, station pressure, wind speed, and wind direction.Source: §3.2.1 · p. 4

3.2.2 STATION ASSESSMENT AND STRATIFICATION

• Station assessment combines temporal alignment and observation-quality checks with source, network, and metadata checks to organize stations into tiers.Released hourly observations are retained within 15 minutes of the target hour; assessment also considers five-variable coverage, seasonal stability, gap lengths, quality-control pass rates, and alignment rates.Source: §3.2.2 · pp. 4–5
• M2Weather organizes 17,173 candidate stations into core, regional augmentation, and extension tiers, then forms Global, Europe, and France settings with 2,504, 983, and 176 stations, respectively.The tiers reflect completeness, stability, metadata consistency, and regional forecasting requirements.Source: §3.2.2 · p. 5

3.2.3 DATA ANALYSIS

• Before imputation, M2-CS median minimum-variable and mean-variable coverage is 96.24% and 98.23%, versus 70.58% and 91.11% for M2-RA; verified identity with current GBON affiliation falls from 35.50% to 24.56% to 19.41% across tiers.The coverage and network evidence together indicate that the tiers capture historical multivariate usability and documented observing-network roles.Source: §3.2.3 · p. 5

3.3 UNIFIED TRAINING AND EVALUATION PROTOCOLS

• Unified protocols use chronological 60%/10%/30% splits, model-specific wrappers that preserve native interaction structures, and common evaluation metrics and test origins.Training uses 48 hours to forecast the next 72 hours; overall metrics are computed in standardized space with equal weight across variables.Source: §3.3 · p. 6

• Backpropagation-free adapters add calibration and enabled interaction corrections to frozen backbone forecasts in dataset-standardized space.The backbone remains fixed in evaluation mode, and the interaction branches process its original predictions in parallel.Source: §4 · p. 6
• Variable interactions mix normalized predictions across variables within each station, using a mapping shared across stations.The mapping does not exchange information between stations.Source: §4 · p. 7
• Station interactions combine one- and two-hop spatial contrasts separately for each variable, without mixing forecast steps.A learnable coefficient sets the strength and direction of each neighborhood correction.Source: §4 · p. 7
• With backbone predictions fixed, squared-error fitting becomes a quadratic problem solved directly by ridge regression for each forecast-step–variable pair.Training-origin–station samples form the feature matrix and residual vector used to fit the coefficients.Source: §4 · p. 7
• The implementation leaves the calibration intercept unpenalized and selects λ by overall standardized validation MSE.It accumulates feature statistics in batches and performs a small regression solve without iterative gradient descent or backpropagation.Source: §4 · p. 7

5 OBSERVATIONS AND FINDINGS

• The evaluation compares 16 representative backbones across four modeling paradigms and three datasets, using 48 hours of history to predict 72 hours.The main text reports mean standardized MSE averaged equally over four target variables.Source: §5 · pp. 7–8, Table 2 · p. 8

5.1 DOES JOINT STATION–VARIABLE MODELING IMPROVE FORECASTING?

• The native MSMV reference, HiSTGNN, has the lowest MSE among 16 backbones across all three datasets, while SSSV has the highest group-mean MSE.On France, group means are 0.4484 for SSSV, 0.3892 for SSMV, 0.3618 for MSSV, and 0.3079 for HiSTGNN; MSSV also beats SSMV at every scale.Source: §5.1 · p. 8, Figure 5 · p. 9

5.2 DOES COMPLETING MISSING INTERACTIONS IMPROVE FORECASTING?

• SSSV, SSMV, and MSSV adapters lower MSE for all 45 model–dataset pairs, with relative reductions of 0.32–19.08%.The group-mean reductions are largest for SSSV at every scale: 10.2% on France, 8.8% on Europe, and 8.4% on Global.Source: Table 3 · p. 8, Figure 5 · p. 9, §5.2 · p. 9
• An interaction-only adapter improves all 30 evaluated model–dataset pairs on France and Global, so gains persist without calibration.For example, Moirai’s France MSE falls from 0.4717 to 0.4094 with interaction-only correction, and to 0.3817 with the full adapter.Source: §5.2 · p. 9

• M2Weather combines quality-controlled datasets at three spatial scales with unified protocols for comparing four station–variable modeling paradigms.The benchmark was introduced to assess the individual and joint roles of station and variable interactions.Source: §6 · p. 9
• Lightweight adapters complement missing interactions in frozen backbone forecasts through closed-form ridge regression without backpropagation.Source: §6 · p. 9
• All 15 adapted backbones achieve lower MSE on all three datasets, and interaction-only ablations retain gains without calibration.The authors identify richer joint architectures, additional forecast horizons, and extreme-weather conditions as directions for extending the analysis.Source: §6 · p. 9

AI USE STATEMENT

• Generative AI assisted preliminary exploration, literature review, and manuscript polishing and editing; the authors retain responsibility for accuracy.Source: AI USE STATEMENT · p. 10

APPENDIX OVERVIEW

• Appendix A covers quality criteria, station tiering and provenance, imputation and target construction, and spatial subsets.Source: APPENDIX OVERVIEW · p. 13
• Appendix B details model hyperparameters and adapter fitting, while Appendix C presents forecasting results, ablations, efficiency analyses, and visualizations.Source: APPENDIX OVERVIEW · p. 13

A.1 QUALITY CRITERIA AND STATION-TIERING THRESHOLDS

• Station records are hourlyized and de-duplicated before coverage and gap statistics; quality-control pass rates and hourly-alignment rates are evaluated separately.Quality metrics include five-year coverage over 43,824 hours, monthly coverage quantiles, gap lengths, weighted QC pass rates, alignment within 15 minutes, and longest station-wide outage.Source: §A.1 · pp. 13–14
• 2,504 M2-CS, 1,637 M2-RA, and 13,032 M2-EP stations result from sequential tiering after quality criteria are applied.Core-tier stations satisfy all M2-CS criteria; after their exclusion, stations satisfying M2-RA criteria form the regional augmentation tier, with the remainder in M2-EP.Source: §A.1 · p. 14

A.2 STATION IDENTITY, PROVENANCE, AND NETWORK MATCHING

• Accepted OSCAR/Surface matches use identifier, location, and supporting metadata evidence; coordinate-only matches within 5 km remain candidates, not official network assignments.The matching procedure also records station identity, location, operating status, and GBON, RBON, and GCOS Surface Network membership.Source: §A.2 · p. 14

A.3 IMPUTATION MASKS AND TARGET-VARIABLE CONSTRUCTION

• Missing intervals of at most 12 hours bounded by valid observations are interpolated; all remaining gaps are filled from the nearest ERA5 grid cell at the same UTC hour.Wind direction is interpolated along the shortest circular arc, correctly handling the 0°/360° boundary; no fallback completion stages were activated in the released datasets.Source: §A.3 · pp. 14–15
• The released data retain five variables, while the forecasting task uses air temperature, derived relative humidity, station pressure, and wind speed.Relative humidity is derived after gap completion from processed air and dew-point temperatures; wind direction remains in the released data, while joint forecasting of both wind quantities is left for future extensions.Source: §A.3 · pp. 15–16

A.4 SPATIAL SUBSET DEFINITIONS AND ADDITIONAL STATISTICS

• Global contains 2,504 M2-CS stations; Europe and France are geographic subsets with 678 and 171 M2-CS stations, respectively.Europe also contains 305 M2-RA stations and France 5; geographic membership does not change a station’s quality tier. Median nearest-neighbor distances are 42.97 km for Global, 32.62 km for Europe, and 24.80 km for France.Source: §A.4 · p. 17

B.1 MODEL HYPERPARAMETERS ACROSS DATASETS

• Trained or fine-tuned backbones share Adam, MSE loss, an initial learning rate of 10^-4, and type1 decay; Timer, TimeMoE, and Moirai use zero-shot evaluation.Training lasts at most 30 epochs, with early-stopping patience 10/5/5 and training-window strides 1/6/6 for France/Europe/Global, so training budgets differ by scale.Source: §B.1 · p. 17
• Model configurations preserve different modeling choices: iTransformer jointly models each station’s four variables, while PatchTST and xPatch process station-variable series independently.Timer and TimeMoE perform zero-shot univariate forecasting; Moirai jointly supplies four variables per station in one inference call.Source: §B.1 · pp. 17–18

B.2 ADAPTER FITTING AND PARAMETER EFFICIENCY

• The station adapter uses a fixed eight-neighbor geographic graph, with two-hop station features computed by applying the same graph twice.The graph uses latitude and longitude only; row-wise softmax of distance-scaled values produces an adjacency matrix that need not be symmetric.Source: §B.2 · p. 19
• Adapter fitting accumulates feature and residual statistics batchwise, then solves one small normal-equation system per forecast step and variable for a fixed λ.The solve avoids iterative gradient updates, backpropagation, and optimizer-state storage; with four target variables, each system has at most seven coefficients.Source: §B.2 · pp. 19–20
• With 72 forecast steps and four variables, station-only, variable-only, and joint adapters contain 1,152, 1,440, and 2,016 coefficients, independent of station count.For iTransformer and TQNet, the 1,152-parameter station adapter is 0.054% and 0.193% of their respective backbone sizes.Source: §B.2 · p. 20, Table 8 · p. 20

C.1 REPORTING CONVENTIONS AND ADAPTER CONFIGURATIONS

• Reported MSE and MAE are per-variable means ± sample standard deviations across seeds 2024, 2025, and 2026; relative changes compare adapted and baseline errors.Negative relative changes indicate improvement, and group-level changes compare equally weighted model means rather than averaging model-wise percentages.Source: §C.1 · pp. 20–21

C.2 SEED-LEVEL SIGNIFICANCE TESTS

• After within-dataset Holm correction, 34 of 36 seed-level comparisons are significant; across datasets, 11/12, 12/12, and 11/12 tests are significant for France, Europe, and Global.The paired tests use three seeds; when all 36 tests form one Holm family, 30 remain significant.Source: §C.2 · p. 21, Table 9 · p. 21

C.3 COMPLETE PER-VARIABLE MSE AND MAE

• All 180 adapted model–dataset–variable combinations improve MSE, while MAE decreases in most but not all displayed comparisons.At displayed precision, MAE decreases in 56/60 France, 54/60 Europe, and 55/60 Global comparisons; the adapter is fitted and selected using squared error.Source: §C.3 · p. 21

C.4 INDIVIDUAL MODELS WITHIN EACH PARADIGM

• Among SSMV models on France, Moirai’s MSE falls 19.08% and Timer-XL’s 4.36%, showing that gains vary by backbone within a paradigm.Figure 6 compares frozen backbones with their full adapters across all 16 backbones.Source: Figure 6 · p. 28, §C.4 · p. 28

C.5 REMOVING THE CALIBRATION BRANCH

• Interaction-only correction lowers MSE for all 15 backbones on France and Global, and the full adapter further improves every evaluated pair.For Global Moirai, MSE falls from 0.4326 to 0.3766 with interaction-only correction and to 0.3585 with the full adapter.Source: §C.5 · p. 28

C.6 SENSITIVITY TO RESIDUAL FITTING SOURCES

• The comparison freezes each backbone and adapter-validation set while varying residual fitting sources across full in-sample, matched in-sample, and subsequent held-out origins.Held-out origins come from a temporal block unseen during backbone training; the original test period remains unchanged.Source: §C.6 · p. 28
• On France, held-out fitting lowers MSE by 7.79%, 7.51%, and 3.19%, with improvements in every seed, and beats equal-sized in-sample matched fitting for all three backbones.In-sample full is slightly better but uses 4.07 times as many origins.Source: Table 23 · p. 30, §C.6 · p. 30

C.7 COMPUTATIONAL EFFICIENCY

• The timing protocol profiles seven backbones on France and Europe at batch size 2, with adapter timing excluding prior backbone training and final testing.Memory is peak GPU allocation for the complete procedure, including the backbone.Source: §C.7 · p. 30
• MSSV models take substantially longer per epoch than SSMV models on both datasets, while native MSMV HiSTGNN is slowest on Europe and second slowest on France.HiSTGNN takes 99.80 minutes on Europe and 118.71 minutes on France; runtime also depends on backbone architecture.Source: Table 24 · p. 30, §C.7 · p. 30
• Complete adapter fitting takes 21.1–45.3% of one backbone training epoch across all 12 cases, including frozen-backbone prediction and validation accumulation.Peak memory is lower than training in 10/12 cases, unchanged for CDPNet on France, and higher for Corrformer on France (28.72 versus 26.31 GiB).Source: Table 25 · p. 31, §C.7 · p. 31

C.8 SENSITIVITY TO LOOK-BACK WINDOW LENGTH

• HiSTGNN has the lowest MSE across all four look-back windows, varying only from 0.3050 to 0.3121, while other backbones’ preferred windows differ.The comparison uses 24-, 48-, 72-, and 96-hour histories with a fixed 72-hour forecast horizon.Source: §C.8 · p. 31, Figure 7 · p. 31

C.9 SENSITIVITY TO KNN NEIGHBORHOOD SIZE

• MSE decreases as K rises from 2 to 16 for all four tested backbones on France and Global, with smaller changes from K = 8 to 16 in most cases.Every configuration improves on its corresponding frozen backbone; the main experiments use K = 8.Source: §C.9 · p. 32, Figure 8 · p. 32

C.10 PREDICTION VISUALIZATIONS

• Adapted forecasts better track warming, humidity decline, wind-speed increase, and gradual pressure variation, though timing, magnitude, rapid peaks, and some fluctuations remain imperfect.The examples cover four Europe forecast windows, each with 48 hours of history and a 72-hour forecast.Source: §C.10 · p. 32, Figure 9 · p. 33, Figure 10 · p. 34