StefaLand: Geoscience Foundation Model

• StefaLand is a geoscience foundation model explicitly designed for dynamic land-surface prediction using global spatiotemporal attributes and forcings.
• It employs a masked autoencoder transformer with cross-variable group masking and residual adapters to improve spatial generalization in data-scarce regions.
• Empirical results show marked improvements over LSTM baselines in streamflow, soil moisture, and soil composition tasks with efficient fine-tuning and physics integration.

Searching arXiv for the specified paper and closely related work to ground the article. arxiv_search(query="(Kraabel et al., 22 Sep 2025)", max_results=5, sort_by="submittedDate")

arxiv_search(query="StefaLand: An Efficient Geoscience Foundation Model That Improves Dynamic Land-Surface Predictions", max_results=10, sort_by="submittedDate")

arxiv_search(query="Geo-Foundation Models: Reality, Gaps and Opportunities", max_results=5, sort_by="relevance")

StefaLand, short for Spatiotemporal Earth Foundation model with Attributes for the Land Surface, is a geoscience foundation model explicitly designed to predict dynamic land-surface interactions—such as streamflow, soil moisture, and soil composition—rather than to recognize patterns in static satellite images. It combines self-supervised pretraining on global spatiotemporal attributes with lightweight, task-specific fine-tuning heads, and is designed to run on academic-level compute. In "StefaLand: An Efficient Geoscience Foundation Model That Improves Dynamic Land-Surface Predictions" (Kraabel et al., 22 Sep 2025), the model is presented as a generative spatiotemporal earth foundation model centered on landscape interactions, with an attribute-based, location-aware architecture and residual fine-tuning adapters intended to improve spatial generalization in data-scarce regions.

1. Problem setting and motivation

Dynamic land-surface variables—streamflow, soil moisture, soil texture, landslides, snow water equivalent, groundwater, vegetation carbon, and related quantities—respond to coupled processes involving heterogeneous static attributes, dynamic forcings, and human influence. The relevant static attributes include topography, soils, geology, land cover, lithology, soil depth, and permeability; the dynamic forcings include precipitation, temperature, radiation, humidity, potential evapotranspiration, and wind; and the human component includes land use change, irrigation, dams, withdrawals, and urbanization. These processes are nonlinear, often thresholded, and can exhibit strong memory, as in soil moisture and groundwater storage (Kraabel et al., 22 Sep 2025).

The model is motivated by three especially challenging issues. The first is spatial generalization and sparsity: high-quality in-situ observations are dense in some developed regions but sparse elsewhere, and machine-learning models trained on dense regions often fail when transferred to data-scarce regions, a setting associated with Prediction in Ungauged Basins (PUB). The second is concept drift across regions: the same input attribute can have different hydrologic meaning depending on climate, soil, and human management, making cross-region parameter transfer unreliable. The third is limited temporal and vertical observability: many crucial variables, including soil moisture at depth, subsurface storage, and soil properties, are not directly observable with satellites at the required resolution and accuracy.

Within this framing, traditional impact models are described as having complementary but limited strengths. Process-based models use mechanistic equations and can generalize physically, but calibration requires local data and parameter regionalization may fail in truly ungauged or conceptually different regions. Statistical and classical machine-learning models are often trained per region or per basin with limited data and use task-specific features. Task-specific deep learning, particularly LSTMs, is strong for hydrologic time series and soil moisture, but is usually trained in a fully supervised fashion on specific datasets and does not exploit the global joint distribution of landscape attributes and forcings.

A central claim of the paper is that recent Earth foundation models built primarily on satellite imagery are ill-suited for dynamic land-surface prediction. The stated reasons are modality mismatch and systems-level mismatch: many key processes and subsurface states are not directly visible; satellite data are noisy, cloud-obscured, and often temporally misaligned with hydrologic processes; and pixel-based imagery does not align with the way hydrologic models typically aggregate information at the basin or site scale. The paper further reports that fine-tuning TerraMind and PrithviWxC on streamflow and soil moisture does not improve over strong LSTM baselines and often underperforms them.

2. Model formulation and architectural design

StefaLand is a masked autoencoder-based transformer that learns to reconstruct masked portions of multivariate spatiotemporal inputs. For a location $s$ and time $t$, the inputs are dynamic forcings $x_{s,t}$ and static landscape attributes $s_s$. During pretraining, the model receives sequences $\{x_{s,t}\}_{t=1}^T$ together with static attributes, randomly masks groups of variables and/or time steps, and reconstructs the masked parts from the unmasked context. The reconstruction target includes both masked dynamic time series and masked static attributes, so the model is generative in the sense of learning to fill in missing parts of the attribute–time–space tensor (Kraabel et al., 22 Sep 2025).

The embedding stage is specified at the level of heterogeneous attributes. For each dynamic variable $c$ at time $t$, the scalar input is embedded by a two-layer MLP with GELU,

$$z_{t,c} = \mathrm{GELU}(x_{t,c} W_{1,c} + b_{1,c}) W_{2,c} + b_{2,c},$$

with $W_{1,c} \in \mathbb{R}^{1 \times 64}$ and $W_{2,c} \in \mathbb{R}^{64 \times 256}$. These embeddings are summed across variables to yield one embedding per time step,

$t$0

Static attributes are embedded similarly and aggregated as a static token, yielding a unified sequence $t$1. A stack of transformer layers then performs multi-head self-attention and feedforward mapping, with positional encodings providing temporal order and the static token enabling every time step to attend to landscape context. Final encoder outputs are linearly projected and passed through a bidirectional LSTM, after which separate decoders reconstruct dynamic and static masked inputs. The pretraining loss is a reconstruction loss over masked slices, normalized per variable and optionally weighted:

$t$2

A distinctive element is Cross-Variable Group Masking (CVGM). Instead of masking individual scalars at random, StefaLand groups variables by physical domain—such as terrain, climate, soil, and vegetation—and masks entire groups jointly across parts of the sequence. This forces the encoder to learn cross-attribute interactions, such as the joint control of runoff and soil moisture by soil texture, slope, and precipitation, or the co-evolution of vegetation indices, land cover, and climate.

The architecture is explicitly location-aware. A static token encodes site-level or basin-level attributes, and temporal tokens attend to it during self-attention. In the notation of the paper, the spatiotemporal input is

$t$3

and the learned temporal and static representations are

$t$4

This design is intended to condition temporal signatures such as flood peaks on elevation, soil, lithology, and related landscape context.

Another defining choice is the use of attribute-based representations rather than raw imagery. Each site or basin is represented by tens of static attributes plus a multivariate forcing time series, which drastically reduces input dimensionality relative to pixel-level image patches. The paper argues that this improves compute efficiency, physical relevance, and global coverage.

For downstream transfer, the pretrained StefaLand encoder is used in frozen form with lightweight adapters and decoders. The principal variant is StefaLand-resConn, in which a task-specific decoder, such as an LSTM, handles temporal sequence prediction while residual paths inject StefaLand features iteratively into the decoder. Other reported variants are StefaLand-no resConn, which uses transformer embeddings once via a simple adapter, and StefaLand-Ablation-resConn, which retains the architecture but uses randomly initialized, non-pretrained StefaLand weights.

3. Pretraining corpus and supervised downstream tasks

The pretraining dataset covers 4,229 locations globally and 40 years of data, roughly 1980–2020. Time-series variables include precipitation, shortwave radiation downwards, relative humidity, max/min temperature, and potential evapotranspiration, drawn from MSWX/MSWEP. Static attributes include land cover fractions and NDVI, soil fractions and soil depth, elevation, slope, aspect, lithology via GLiM, hydrogeology via GLHYMPS, socio-economic variables such as population density and GDP, and forest intactness. Static and dynamic datasets are combined into unified samples per location, with one sample defined as a site-level static attribute vector together with a 40-year forcing time series. Data are stored in shards and lazily loaded, and shards are randomly reindexed at each epoch to ensure diverse sampling without placing all data in memory (Kraabel et al., 22 Sep 2025).

Downstream evaluation is organized around three main tasks and four datasets. The first task is streamflow prediction on CAMELS (US), using 531 basins with clear watershed boundaries. Two generalization regimes are used: PUB, implemented as 10 random folds with leave-one-fold-out, and PUR, implemented as 7 contiguous regions with regional leave-one-out. The compared models include LSTM-SL, StefaLand-Transformer trained from scratch on task data, StefaLand-resConn, StefaLand-no resConn, StefaLand-Ablation-resConn, TerraMind-resConn, and physics-based variants such as LSTM-HBV1.1 and StefaLand-resConn-HBV1.1.

The second streamflow benchmark is global headwater basins, comprising 3,434 basins filtered for data completeness. Experiments include a random holdout setting and a regional holdout across continents (RH-C) in which all but one continent are used for training and the excluded continent is used for testing. The paper reports results for North America, South America, and Europe in this cross-continent setup.

The third task is soil moisture prediction using the International Soil Moisture Network (ISMN), with 1,316 sites globally. Dynamic inputs include precipitation from GPM and MSWEP, ERA5-Land forcings, albedo from MODIS, and land surface temperature. Static attributes include elevation, slope, aspect, roughness, curvature, soil properties from HWSD, land cover, and NDVI. The evaluation uses random location holdout with 5-fold spatial cross-validation and a regional holdout in which all continents except Europe are used for training and 129 European sites are used for testing. Models include LSTM-SL, StefaLand-Transformer, StefaLand-resConn, StefaLand-no resConn, Ablation-resConn, TerraMind-resConn, and PrithviWxC-resConn.

The fourth dataset supports soil composition estimation from ISRIC WoSIS profiles. The goal is to predict in-situ soil texture fractions—clay, sand, and silt—using StefaLand pretrained on HWSD soil attributes. For each property, complementary fractions are masked to avoid leakage; inputs include climate, terrain, land cover, and other soil/environment variables; and the compared models are StefaLand-Transformer, StefaLand-Finetune, and a random forest baseline.

4. Empirical performance and transfer behavior

The reported empirical results are centered on spatial generalization. On CAMELS PUB, LSTM-SL achieves RMSE 1.402, Corr 0.762, NSE 0.636, while StefaLand-resConn achieves RMSE 1.111, Corr 0.869, NSE 0.717, corresponding to an approximately 20% reduction in RMSE relative to the LSTM baseline. StefaLand-no resConn attains RMSE 1.171, Corr 0.823, NSE 0.706. TerraMind-resConn records RMSE 1.332, Corr 0.777, NSE 0.637, described as showing no real gain over the LSTM baseline. Ablation-resConn yields RMSE 1.355, Corr 0.801, NSE 0.661, which the paper interprets as evidence that architecture alone is not sufficient and that pretraining matters (Kraabel et al., 22 Sep 2025).

On CAMELS PUR, LSTM-SL obtains RMSE 1.609, Corr 0.743, NSE 0.554, whereas StefaLand-resConn obtains RMSE 1.344, Corr 0.801, NSE 0.635, corresponding to an approximately 17% RMSE reduction and a substantial improvement in NSE. The paper further states that the physics-based StefaLand-resConn-HBV1.1 improves NSE in PUR, indicating a synergy between physics and pretrained features.

On the global streamflow benchmark, StefaLand is again reported to outperform supervised LSTM baselines. In random holdout, LSTM-SL achieves RMSE 0.870, Corr 0.798, while StefaLand-Transformer achieves RMSE 0.749, Corr 0.843, or approximately 14% lower RMSE with higher correlation. TerraMind performs substantially worse at RMSE 1.156, Corr 0.580. In RH-C, LSTM-SL records RMSE 1.253, Corr 0.672, and StefaLand-Transformer records RMSE 1.075, Corr 0.697, corresponding to an approximately 14.1% RMSE reduction. The paper notes that µbRMSE is close to RMSE for StefaLand, suggesting that errors are mostly random rather than systematic bias.

For soil moisture on ISMN, the strongest reported model is again StefaLand-resConn. In random holdout, LSTM-SL reaches RMSE 0.073, Corr 0.764, while StefaLand-resConn reaches RMSE 0.068, Corr 0.783. StefaLand-no resConn records RMSE 0.075, Corr 0.741. TerraMind-resConn and PrithviWxC-resConn are weaker, at RMSE 0.083, Corr 0.694 and RMSE 0.081, Corr 0.703, respectively. In the Europe regional holdout, LSTM-SL yields RMSE 0.112, Corr 0.510, whereas StefaLand-resConn yields RMSE 0.090, Corr 0.638, described as an approximately 25% improvement in correlation over the LSTM baseline.

For soil composition estimation, the reported pattern differs by target variable. StefaLand-Finetune improves prediction of clay and sand relative to the random forest baseline: for clay, it reaches Corr 0.51, $t$5 versus Corr 0.05, $t$6 for random forest; for sand, it reaches Corr 0.70, $t$7 versus Corr 0.26, $t$8 for random forest. For silt, however, random forest performs better at Corr 0.84, $t$9, compared with Corr 0.63, $x_{s,t}$0 for StefaLand-Finetune. The paper suggests that this may reflect the way silt is calculated, often as the remainder of sand plus clay, together with dataset noise.

Task and setup	Baseline result	StefaLand result
CAMELS PUB streamflow	LSTM-SL: RMSE 1.402, Corr 0.762, NSE 0.636	StefaLand-resConn: RMSE 1.111, Corr 0.869, NSE 0.717
CAMELS PUR streamflow	LSTM-SL: RMSE 1.609, Corr 0.743, NSE 0.554	StefaLand-resConn: RMSE 1.344, Corr 0.801, NSE 0.635
Global RH-C streamflow	LSTM-SL: RMSE 1.253, Corr 0.672	StefaLand-Transformer: RMSE 1.075, Corr 0.697
ISMN Europe soil moisture	LSTM-SL: RMSE 0.112, Corr 0.510	StefaLand-resConn: RMSE 0.090, Corr 0.638

These results support the paper’s central argument that pretraining on globally distributed attributes and forcings can improve PUB, PUR, and cross-continent transfer relative to strong supervised baselines. They also support the narrower architectural claims that residual integration is consistently stronger than one-shot adapters and that image-based foundation models are not automatically effective when repurposed for dynamic land-surface prediction.

5. Generalization mechanisms and computational profile

The paper attributes StefaLand’s spatial generalization to four mechanisms. The first is foundation model pretraining on 40 years of globally distributed attributes, which is said to encourage physically relevant patterns rather than dataset-specific quirks. The second is CVGM, which forces learning of physical relationships across climate, terrain, soil, and vegetation domains. The third is the use of residual adapters with a frozen encoder, which is argued to preserve global knowledge and reduce overfitting to data-rich regions. The fourth is physics-based integration, especially via HBV1.1, which adds process constraints and further improves generalization in some settings (Kraabel et al., 22 Sep 2025).

The computational profile is presented as a major design feature. The transformer backbone has approximately 12 million parameters. Pretraining requires approximately 16,000 seconds per epoch on 6× NVIDIA V100 GPUs (240 GB total) and about 720 GPU hours in total. Fine-tuning StefaLand-resConn requires approximately 30 seconds per epoch on 2× V100 (80 GB), compared with approximately 12 seconds per epoch for the LSTM baseline. Physics-based HBV variants require 280–320 seconds per epoch, reflecting the additional expense of physics integration.

These figures are contrasted with several image-centric Earth foundation models: TerraMind at approximately 7,680 GPU hours, Aurora at approximately 14,592 GPU hours, and AlphaEarth at approximately 28,672 GPU hours. The paper attributes StefaLand’s efficiency to three factors: attribute-based representation, masked-autoencoder self-supervision, and the use of a frozen encoder with small downstream decoders. A plausible implication is that the model’s accessibility is not incidental but built into the methodological design, since academic groups with modest GPU clusters are explicitly identified as intended users.

6. Applications, limitations, and nomenclature

StefaLand is presented as a general-purpose land-surface foundation layer. The stated application areas include hydrology, where it can support streamflow forecasting, regional runoff mapping in ungauged basins and continents, and hybrid physics-informed models such as HBV1.1 parameterization; soil and agriculture, where it can support soil moisture monitoring and forecasting, soil property inference and harmonization across datasets such as HWSD, ISRIC, SoilGrids, and national soil maps, as well as estimation of root-zone storage and plant-available water capacity; geohazards, including landslide susceptibility, wildfire risk, and future extension to snow water equivalent; and climate impacts and resource management, including climate-driven changes in hydrologic regimes, water allocation, infrastructure planning, and hazard mitigation in data-scarce regions (Kraabel et al., 22 Sep 2025).

The paper also states that code and pretrained weights are publicly released and that the implementation includes full pretraining and fine-tuning pipelines together with detailed variable lists. This suggests an intended role not only as a benchmark model but also as a reusable encoder in differentiable workflows, including integration with process-based models and domain-specific task heads.

Several limitations are explicit. Pretraining is focused on selected variables, mostly surface climate, soils, terrain, and land cover; deep groundwater and complex human systems are not directly represented. Benchmarking still relies heavily on regions with relatively good data availability, especially the US and Europe, so performance in extremely data-scarce areas such as Central Africa and remote parts of Asia remains to be fully validated. The model is currently attribute-based, and future work is described in terms of adding higher-frequency forcings, anthropogenic and ecosystem dynamics such as urban growth, deforestation, and irrigation, image-like inputs for tasks where high-resolution spatial patterns matter, sequence-aware pretraining objectives, uncertainty-aware heads, and improved integration with differentiable physics.

A separate nomenclature issue is relevant. In an exposition of "Approximation of the interface condition for stochastic Stefan-type problems" (Mueller, 2018), the term “StefaLand” is used informally to denote the “Stefan-problem landscape.” That usage concerns stochastic moving boundary problems, semilinear SPDEs, and approximation of Stefan-type interface conditions by local volume imbalance. It is distinct from the land-surface foundation model discussed here. The coexistence of these usages can create terminological ambiguity, but the geoscience literature identified by (Kraabel et al., 22 Sep 2025) uses StefaLand specifically for the attribute-based Earth foundation model for dynamic land-surface prediction.