A Physics-Guided Spatio-Temporal Bidirectional LSTM Model for Systemic Nighttime LST Forecasting and Trend Projection: A Case Study of Afghanistan

1. Introduction

Land surface temperature (LST) is a fundamental geophysical variable that regulates land–atmosphere interactions, surface energy balance, hydrological processes, and ecosystem dynamics [1]. In the context of accelerating global warming, accurate monitoring and forecasting of LST have attracted growing attention from both the scientific community and operational climate-management agencies [2]. Within this body of research, daytime LST has been shown to respond strongly to instantaneous solar radiation and short-term surface conditions, whereas nighttime LST (LSTn) more directly reflects land-surface thermal inertia, subsurface heat storage, and background atmospheric stability [3]. Reliable LSTn forecasting is therefore particularly valuable for climate assessment, environmental risk evaluation, and the design of adaptation strategies in climatically vulnerable regions [4]. A wide range of traditional statistical and deep learning techniques has been employed in the related literature to identify the dominant drivers of LSTn and to forecast its spatiotemporal fluctuations. Traditional approaches, including Multiple Linear Regression (MLR), Geographically Weighted Regression (GWR), Autoregressive Integrated Moving Average (ARIMA), Seasonal ARIMA (SARIMA), Support Vector Regression (SVR), Random Forest (RF), Extreme Gradient Boosting (XGBoost), and the Mann-Kendall trend test coupled with Sen’s slope estimator—have been widely applied to characterize long-term LSTn trends and to quantify the influence of biophysical and climatic covariates [5,6]. More recently, deep learning (DL) architectures have demonstrated superior capability in capturing the nonlinear and spatiotemporal dependencies inherent in LST dynamics. Notable examples include Artificial Neural Networks (ANN), Long Short-Term Memory (LSTM) networks, Bidirectional LSTM (BiLSTM), Gated Recurrent Units (GRU), Convolutional Neural Networks (CNN), ConvLSTM, Transformer-based models, and hybrid frameworks such as CNN-LSTM and Transformer-LSTM, which integrate spatial feature extraction with temporal sequence modeling to enhance forecasting accuracy [7,9]. Among these DL techniques, Physics-guided Spatio-temporal Bidirectional LSTM (PG-ST-BiLSTM) has emerged as a promising approach in which domain knowledge is embedded into data-driven architectures to improve predictive consistency [10,11]. Therefore, this study employs PG-ST-BiLSTM to systematically evaluate the accuracy of the selected method on historical data spanning from 2001 to 2024, and subsequently to forecast LSTn from 2025 to 2040 using monthly data from Afghanistan.

To realize this forecasting framework, the methodological workflow typically encompasses several interdependent stages, including the acquisition of LSTn and environmental driver datasets, time-series preprocessing and feature construction, development and training of the forecasting model, generation of temporal predictions, and subsequent validation [12,13]. Considerable effort has been devoted to the construction and optimization of statistical and DL models for LSTn estimation and forecasting. For instance, LSTM- and Transformer-based architectures have been refined through hyperparameter tuning and cross-validation strategies to enhance temporal accuracy [14,15], while advanced optimization techniques such as Bayesian optimization, adaptive learning-rate scheduling, and regularization have been adopted to improve model stability and generalization in environmental time-series applications [16,17].

However, these model-centric improvements are only one side of the forecasting pipeline. Equally critical-yet far less systematically addressed-is the design of the temporal input dataset itself. In many cases, the length of the input sequence, the representation of seasonal cycles, and the integration of multivariate environmental drivers are selected empirically rather than through principled evaluation. Short input windows may fail to capture full seasonal periodicity, whereas excessively long sequences may introduce noise and unstable gradient propagation [18,19]. Similarly, the temporal alignment and joint representation of LSTn with vegetation, precipitation, and soil moisture drivers are often implemented without explicit consideration of cross-seasonal thermal persistence. As a result, models tend to learn short-term correlations while underrepresenting the multi-month memory effects that govern nighttime surface thermal dynamics. Moreover, most forecasting frameworks rely on pointwise loss functions that optimize prediction error independently at each time step, without constraining the continuity of month-to-month temperature transitions [20].

Despite these advances, existing studies have predominantly emphasized architectural refinement and parameter optimization [21,22], yet the systematic design of the temporal input dataset including feature selection, gap-filling, and temporal alignment-has received limited attention, especially in data-sparse regions such as Afghanistan. This limitation is particularly problematic in regions characterized by strong elevation gradients and nonlinear seasonal variability, where physically inconsistent forecasts can readily emerge due to the neglect of terrain-driven thermal dynamics. Consequently, despite ongoing improvements in model architecture and hyperparameter optimization, achieving temporally coherent and physically consistent nighttime LSTn forecasting across complex mountainous environments remains a significant challenge.

Many of the limitations identified above stem from methodological choices that are made empirically rather than through systematic evaluation, introducing avoidable uncertainty into LSTₙ predictions [5,8]. Although uncertainty in satellite-based LSTn forecasting has attracted increasing attention in recent years, existing studies have typically examined individual sources of variability in isolation. Preprocessing decisions-including quality-control filtering, temporal aggregation (e.g., 8-day to monthly composites), and gap-filling of cloud- or snow-flagged Moderate Resolution Imaging Spectroradiometer (MODIS), observations have been shown to alter the multi-month thermal signal that the model is required to learn [5,6,8,23], while spatial resolution and resampling strategies further affect the representation of thermal gradients in heterogeneous mountainous regions [24,25]. Equally important, the configuration of temporal inputs and environmental drivers, including input-window length, seasonal representation, and multivariate driver integration, is often determined through prior experience or trial-and-error, even though spatial heterogeneity and strong elevation, climate gradients, such as those characterizing Afghanistan, mean that a configuration well suited to one basin may be suboptimal in another [18,19,26]. At the architectural level, recurrent and convolutional networks cannot, on their own, prioritize the distant temporal dependencies that govern cross-seasonal nocturnal thermal memory; Transformer-based attention has been introduced to overcome this limitation [27], and BiLSTM architectures have been shown to improve bidirectional contextual learning across seasonal cycles [28,29]. Nevertheless, no prior study has systematically benchmarked a comprehensive set of DL architectures, including Temporal Convolutional Networks (TCN), Long Short-Term Memory (LSTM), and Bidirectional LSTM (BiLSTM), against a PG-ST-BiLSTM within a unified LSTn forecasting framework.

By effectively addressing these limitations, this study makes several key contributions to the literature. First, this study proposes the PG-ST-BiLSTM framework-a unified pipeline that integrates temporal input design, a hybrid Transformer-BiLSTM architecture, and a physics-guided training objective to improve both the predictive accuracy and physical consistency of LSTn forecasting in data-scarce mountainous environments. Second, three environmental drivers are employed as predictors, including NDVI, precipitation (PREC), and soil moisture (SM), selected for their collective capacity to capture vegetation phenology, atmospheric water input, and subsurface thermal inertia, which are the dominant processes governing nocturnal surface energy exchange in Afghanistan’s continental climate. Third, rather than optimizing a single model in isolation, this study systematically evaluates three methodological factors that are typically treated independently in the literature: (1) the configuration of the temporal input window and multivariate driver representation; (2) the integration of Transformer-based long-range attention with bidirectional recurrent seasonal learning; and (3) the incorporation of a physics-motivated temporal smoothness constraint directly into the loss function. Each factor is evaluated against four baselines, including TCN, LSTM, BiLSTM, and Linear Regression, under identical data partitions and preprocessing conditions, enabling performance differences to be attributed to specific design choices rather than to data or implementation variation. Finally, the novelty of PG-ST-BiLSTM lies in its principled integration of components, which can be articulated through three distinct contributions that substantively differentiate this study from prior work: (1) the first controlled benchmarking of TCN, LSTM, BiLSTM, and PG-ST-BiLSTM on LSTn, providing direct empirical evidence of which architectural choices matter most; (2) a physics-motivated smoothness constraint derived from land-surface thermal inertia principles, demonstrated through ablation to contribute more to predictive accuracy than either the Transformer encoder or the multi-driver input set; and (3) the first basin-disaggregated, season-stratified nighttime warming projections through 2040 for Afghanistan, a region characterized by extreme topographic and climatic heterogeneity, near-zero ground-based monitoring infrastructure, and high practical stakes for water resource planning, agricultural risk assessment, and cryosphere monitoring.

2. Study area and data sources

2.1. Study area

This study performs spatially heterogeneous LSTn forecasting across the five major river basins of Afghanistan, including Amu Darya, Harirud–Murghab, Helmand, Kabul, and Northern. Afghanistan is located in south-central Asia between 29°-38° N and 61°-74° E (Figure 1) [32]. The country covers approximately 652,000 km² and exhibits substantial topographic variability, with elevations ranging from approximately 230 m in the southwestern lowlands to 7,471 m in the Hindu Kush Mountains in the northeast. This pronounced elevational gradient generates strong climatic contrasts across relatively short horizontal distances: northeastern high-altitude regions experience cold winters and comparatively high annual precipitation, whereas the southern lowlands are characterized by hot, arid conditions with limited rainfall. Hydrologically, the country is drained by five aforementioned major river basins, each exhibiting distinct hydroclimatic regimes [33,34]. The combined spatial heterogeneity in elevation, climate, and land cover gives rise to complex nocturnal surface thermal dynamics, making Afghanistan a particularly informative setting in which to investigate LSTn variability for regional climate assessment and environmental management [35,36].

2.2. Nighttime Land Surface Temperature Dataset

The LSTn data constitute the primary variable analysed in this study. LSTn was derived from the MODIS MOD11A2 Collection 6.1 product acquired from the Terra satellite platform [37]. This dataset provides 8-day composite LSTn observations at a spatial resolution of 1 km for the period 2001–2024 [38]. The LST_Night_1km band was selected to represent nocturnal surface thermal conditions. MODIS LST digital numbers were converted to degrees Celsius using the standard scale factor: $LST(°C)=DN\times 0.02-273.15\left[39\right]$. Quality control filtering was applied using the QC_Night band to remove cloud-contaminated and low-quality pixels. Pixels flagged as good quality under the mandatory QA bits and acceptable LST error conditions were retained for analysis. A total of 1,103 valid 8-day composite images were processed for the study period. The mean and trend of LSTn, NDVI, SM, and PREC are visualized in Figure 2(a)-2(h), respectively.

2.3. Environmental Driver Datasets

Based on previous studies of land–atmosphere interactions and surface energy balance processes, and considering data availability and consistency over the study period, three environmental driver datasets were incorporated into the analysis. An overview of the datasets is provided in Table 1.

2.4. Vegetation Index (NDVI)

Vegetation dynamics were represented using the MODIS MOD13Q1 Collection 6.1 Normalized Difference Vegetation Index (NDVI) product. This dataset provides 16-day composite observations at 250 m spatial resolution for 2001–2024 [39]. NDVI captures vegetation greenness and canopy activity, which influence surface energy partitioning and nocturnal cooling behaviour [41]. The NDVI data were resampled to 1 km resolution to match the MODIS LST grid.

2.5. Precipitation

Precipitation data were obtained from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) dataset. CHIRPS provides daily rainfall estimates at approximately 0.05° (~5.5 km) spatial resolution for the study period [42]. Daily precipitation values were aggregated into monthly totals and subsequently resampled to the 1-km LST grid. Precipitation is a key climatic driver that influences soil moisture conditions, vegetation growth, and surface thermal properties [43,44].

ECMWF Reanalysis v5 (ERA5) reanalysis was considered as an alternative source for both precipitation and soil moisture. However, CHIRPS was preferred over ERA5 precipitation for three reasons: (1) CHIRPS merges satellite-derived infrared estimates with quality-controlled ground station data, making it more accurate than pure reanalysis in data-sparse regions such as Afghanistan where ERA5’s precipitation fields are poorly constrained by in-situ observations [45]; (2) CHIRPS provides a higher native spatial resolution (~5.5 km vs ERA5’s ~31 km) [46], which better captures the orographic precipitation gradients characteristic of Afghanistan’s mountainous terrain; and (3) CHIRPS has been independently validated as the most reliable gridded precipitation product for Central Asia and the Hindu Kush region in multiple benchmark studies [47,48]. Similarly, FLDAS NOAH01 was preferred over ERA5 soil moisture because FLDAS assimilates CHIRPS precipitation as its forcing input, ensuring internal consistency between the precipitation and soil moisture fields used in this study, a consistency that would be compromised if ERA5 soil moisture were paired with CHIRPS precipitation.

2.6. Soil Moisture (0–10 cm)

Near-surface soil moisture data were derived from the FLDAS (NOAH01 V001) dataset. This product provides monthly soil moisture estimates for the 0–10 cm soil layer at approximately 0.1° (~11 km) spatial resolution for 2001-2024 [49]. Soil moisture strongly affects land surface thermal inertia and nocturnal heat storage capacity [50]. The data were spatially resampled to the 1-km LST grid to ensure consistency with other variables.

3. Methodology

3.1. PG-ST-BiLSTM Framework and Key Steps for LSTn Modeling

The common nighttime LST forecasting framework generally includes the following processes (Figure 3): (1) obtain satellite-derived LST data and associated environmental driver datasets; (2) conduct systematic data preprocessing and feature construction; (3) transform the time series into supervised learning sequences; (4) select appropriate forecasting models to construct the prediction framework; (5) perform short- to medium-range forecasting; and (6) evaluate prediction accuracy and quantify forecast uncertainty using statistical metrics. Each of these steps may introduce methodological uncertainty if a single modeling configuration is subjectively selected.

Therefore, this study proposes a structured Physics-Guided Spatio-Temporal LSTM (PG-ST-LSTM) framework to ensure methodological robustness and physical consistency (Figure 4) [9,51]. First, LSTn and associated environmental drivers were obtained from MODIS MOD11A2 Collection 6.1 (8-day, 1 km; 2001-2024). Three predictor variables were integrated: NDVI (MOD13Q1), precipitation (CHIRPS), and near-surface soil moisture (FLDAS NOAH01). These variables represent vegetation dynamics, hydrological forcing, and surface moisture conditions influencing nocturnal thermal behaviours [42,43,49,50]. Second, data preprocessing was performed to ensure spatial and temporal consistency. MODIS QA flags were applied to remove cloud-contaminated pixels [52]. The 8-day LST composites were aggregated to monthly means, and all ancillary datasets were resampled to the 1 km grid. Basin-level zonal statistics and seasonal (DJF, MAM, JJA, SON) averages were computed. All predictors were normalized using Min-Max scaling to stabilize model training. Third, the multivariate time series was converted into supervised learning sequences using a 12-month sliding window to capture the annual thermal cycle.

Each input vector comprised six features: LSTn, NDVI, precipitation, soil moisture, and cyclic seasonal encodings (sin, cos). For the validated short-to-medium horizon (2015-2024), this direct multi-output strategy avoids recursive error accumulation by predicting all four steps simultaneously. For the long-range projection period (2025-2040), an autoregressive rolling strategy was employed, wherein the model’s four-step output is appended to the input sequence and the window is advanced by one step iteratively.

This transition is deliberately conservative: the direct strategy is used where ground-truth comparison is possible, while the autoregressive strategy is adopted only for the projection horizon where no observed data exist for verification, with forecast uncertainty explicitly quantified through the confidence-tier framework described in Section 3.4, and chronological partitioning (2001-2014 training; 2015-2024 validation) was applied to prevent temporal leakage. Fourth, suitable forecasting models were implemented to construct the prediction framework. In addition to four baselines (LR, LSTM, BiLSTM, TCN), the proposed PG-ST-LSTM model was developed.

The architecture integrates a dual-block Transformer encoder for long-range temporal attention, stacked Bidirectional LSTM layers for sequential learning, and a physics-guided temporal smoothness constraint within the loss function to enforce physically plausible trajectories. The ‘spatio-temporal’ designation in PG-ST-BiLSTM refers to the framework’s simultaneous operation across five spatially distinct river basins under a unified modeling protocol, wherein basin identity is encoded implicitly through basin-specific input sequences rather than through explicit spatial convolution or graph-based spatial attention.

This is consistent with the use of ‘spatio-temporal’ in multi-site time series forecasting studies [53,54], where the term denotes joint modeling across spatial units rather than pixel-level spatial feature extraction. Readers seeking explicit spatial convolution should note that the basin-level aggregation performed during preprocessing converts the 2D spatial domain into a set of representative time series, making pixel-level CNN-based spatial processing neither necessary nor appropriate for this forecasting task.

Fifth, short- to medium-range forecasting was conducted using direct multi-output prediction to avoid recursive error accumulation. The model was trained using the Adam optimizer with early stopping and learning rate decay. Validation was performed on the independent 2015-2024 period to assess temporal generalization.

Sixth, prediction accuracy and uncertainty were quantified using R², RMSE, and MAE. Seasonal performance was examined to evaluate robustness under varying climatic conditions. Ablation experiments were conducted to isolate the contributions of the Transformer encoder, environmental drivers, and physics-guided regularization. Forecast degradation was assessed empirically (0.069 °C/year), forming the basis for confidence-tier classification in long-range projections. Mann–Kendall and Sen’s slope analyses were subsequently applied to derive statistically robust trend-constrained projections beyond the validated forecasting horizon [55,56].

3.2. Dataset Preparation for LSTn Modeling

3.2.1. Multicollinearity Assessment

The assessment of multicollinearity among environmental predictors is an important step in LST forecasting, particularly for regression-based baselines and correlation analysis. Multicollinearity arises when predictor variables exhibit strong interdependence, such that variation in one variable is systematically associated with variation in another. This redundancy may lead to unstable coefficient estimation, reduced interpretability, and inflated variance in linear modeling frameworks.

To ensure the statistical independence of the input predictors and to avoid redundancy in the DL model, the Variance Inflation Factor (VIF) was computed for each LSTn predictor following [57,58]. The VIF is defined as:

$$VI{F}_i={\displaystyle \frac{1}{1-R_i^2}}$$

(1)

where $R_i^2$ represents the coefficient of determination obtained by regressing the predictor $i$ against all remaining predictors. The reciprocal term $1-R_i^2$, referred to as tolerance, provides an equivalent measure of collinearity. A VIF value exceeding commonly accepted thresholds (e.g., VIF > 5 or 10) indicates potential multicollinearity [59]. Only predictors satisfying acceptable VIF criteria were retained for regression-based baseline modeling and driver analysis. Although DL architectures such as PG-ST-LSTM are generally less sensitive to multicollinearity due to nonlinear representation learning, this assessment ensures statistical robustness and interpretability in baseline comparisons and correlation analyses.

3.2.2. Temporal Feature Construction for Continuous Predictor Variables

Feature representation in LST forecasting requires careful design to capture the nonlinear relationship between environmental drivers and nighttime temperature. Unlike traditional regression-based discretization approaches, this study employs continuous predictor variables processed through engineered temporal encodings. Rather than dividing predictors into discrete intervals, seasonal periodicity is explicitly encoded using cyclical sine and cosine transformations.

The input feature representation consists of six continuous variables:

➢

LSTn (target variable);

➢

NDVI (vegetation state);

➢

Precipitation (hydrological forcing)

➢

Soil moisture (surface moisture availability);

➢

Sine-transformed month (seasonal phase);

➢

Cosine-transformed month (seasonal phase).

The cyclical encoding is formulated as:

$${\mathrm{Month}}_{sin}=sin\left({\displaystyle \frac{2\pi \times \mathrm{month}}{12}}\right),{\mathrm{Month}}_{cos}=cos\left({\displaystyle \frac{2\pi \times \mathrm{month}}{12}}\right)$$

(2)

This approach preserves annual thermal continuity, ensures December and January are treated as adjacent, and eliminates artificial discontinuities inherent in discrete interval methods. All predictors were normalized using Min-Max scaling to the [0,1] range, standardizing gradient magnitudes, and ensuring numerical stability during training [60,61].

3.2.3. Pearson Correlation and Driver Association Analysis

To quantify the biophysical influences on nighttime LST, Pearson correlation coefficients were computed between monthly basin-mean LSTn and the three environmental drivers (NDVI, precipitation, and soil moisture). The Pearson correlation coefficient measures the linear association strength between two variables and is formulated as [62]:

$$r={\displaystyle \frac{{\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}$$

(3)

Correlations were assessed both annually and by season (DJF, MAM, JJA, SON) to identify whether driver–LSTn relationships vary across the annual cycle. A correlation coefficient magnitude exceeding a threshold of r > 0.7 indicates a strong association, while r < 0.3 suggests weak influence. Analyses were conducted over the complete 2001-2024 time series (n = 288 months) and tested for statistical significance at p < 0.05. Following Pearson analysis, multiple linear regression was applied to quantify the combined explanatory power of the three drivers:

$$LSTn={\beta }_0+{\beta }_1\cdot NDVI+{\beta }_2\cdot Precip+{\beta }_3\cdot SoilM+ϵ$$

(4)

The resulting seasonal R² values serve as a linear predictability baseline against which PG-ST-LSTM nonlinear performance is compared. This enables quantification of the performance gain achievable through deep learning architectures.

3.2.4. Sequence Construction and Chronological Sample Partitioning

To prevent information leakage from future observations into model training, strict chronological partitioning was applied. The complete 2001-2024 time series (288 months) was divided as follows:

Training period: January 2001- December 2014 (156 input-output sequence pairs)

Validation period: January 2015 - December 2024 (116 sequence pairs)

A sliding-window approach with a step size of one month was employed to construct sequinces, such that consecutive sequences overlap by 11 input months. Each sequence consists of a 12-month input window and a 4-month direct output vector. This 12-month input window provides the model with one complete annual thermal cycle, the dominant mode of LSTn variability in Afghanistan’s continental climate. The dataset imbalance was avoided by maintaining consistent sample representation across seasons. Both training and validation periods span multiple complete years, ensuring representative seasonal coverage in both subsets.

3.3. Machine Learning Models for Nighttime LST Forecasting

3.3.1. Linear Regression

Linear Regression (LR) serves as a fundamental baseline in this study. LR models the relationship between predictors and the target variable through a linear function [63]:

$$\widehat{LSTn}={\beta }_0+{\displaystyle \sum _{i=1}^p}{\beta }_i{X}_i$$

(5)

where ${\beta }_i$ are regression coefficients and $X_i$ are predictor variables. Despite its simplicity, LR provides interpretable coefficient estimates and establishes a minimal performance benchmark. Limitations include its inability to capture nonlinear temporal dependencies and interactions among environmental drivers. LR serves as the lower-bound baseline for assessing deep learning performance gains.

3.3.2. Long Short-Term Memory (LSTM)

The Long Short-Term Memory (LSTM) [64] It is a recurrent neural network architecture designed to capture temporal dependencies in sequential data. Unlike standard feedforward networks, LSTM maintains hidden states across time steps and employs gating mechanisms [65] to regulate information flow:

$$f_t=\sigma (W_f\cdot [h_{t-1},x_t]+b_f)$$

(6)

$$i_t=\sigma (W_i\cdot [h_{t-1},x_t]+b_i)$$

(7)

$${\tilde{c}}_t=tanh(W_c\cdot [h_{t-1},x_t]+b_c)$$

(8)

$$c_t=f_t\odot c_{t-1}+i_t\odot {\tilde{c}}_t$$

(9)

$$o_t=\sigma (W_o\cdot [h_{t-1},x_t]+b_o)$$

(10)

$$h_t=o_t\odot tanh\left(c_t\right)$$

(11)

Te forget gate $f_t$, input gate $i_t$, and output gate $o_t$ control the flow of information. LSTM effectively captures long-range temporal patterns and has demonstrated success in time-series forecasting. However, standard unidirectional LSTM lacks explicit mechanisms for attending to long-range temporal dependencies across the input window.

3.3.3. Bidirectional LSTM (BiLSTM)

Bidirectional LSTM [29,66] processes sequences in both forward and backward directions, enabling the model to exploit contextual information from both past and future observations within the input window:

$${\overrightarrow{h}}_t={\mathrm{LSTM}}_{\mathrm{forward}}(x_t,{\overrightarrow{h}}_{t-1})$$

(12)

$${\overleftarrow{h}}_t={\mathrm{LSTM}}_{\mathrm{backward}}(x_t,{\overleftarrow{h}}_{t+1})$$

(13)

$$h_t=[{\overrightarrow{h}}_t;{\overleftarrow{h}}_t]$$

(14)

This bidirectional context is particularly valuable for capturing asymmetric warming and cooling transitions in continental climates. BiLSTM improves upon unidirectional LSTM by leveraging information from the complete input sequence. However, it still lacks explicit long-range attention mechanisms for prioritizing distant temporal relationships.

3.3.4. Temporal Convolutional Network (TCN)

The Temporal Convolutional Network (TCN) [67,68] applies dilated causal convolutions to extract temporal patterns from sequential data. Unlike recurrent architectures, TCN processes entire sequences through parallel convolutional operations, enabling efficient parallel computation:

$$y_t=f(x_t,x_{t-1},\dots ,x_{t-R+1})$$

(15)

where $R$ The receptive field is determined by dilation factors and kernel size. TCN can capture long-range dependencies through exponentially increasing dilation rates without the gradient vanishing issues of RNNs. However, TCN lacks explicit attention mechanisms and may not efficiently capture multi-scale temporal interactions present in seasonal climate data.

3.3.5. Physics-Guided Spatio-Temporal LSTM (PG-ST-LSTM)

The proposed PG-ST-LSTM integrates three architectural innovations to overcome identified limitations:

Transformer Encoder: A dual-layer Transformer encoder [7] With multi-head self-attention (8 heads, 64-dimensional keys) captures long-range dependencies regardless of temporal distance:

$$\mathrm{Attention}(Q,K,V)=\mathrm{softmax}\left({\displaystyle \frac{Q{K}^T}{\sqrt{d_k}}}\right)V$$

(16)

Bidirectional LSTM Layers: Two stacked BiLSTM layers [29,66] (64 and 32 units) process Transformer-encoded sequences in both directions, exploiting bidirectional seasonal context.

Physics-Guided Loss Function: A physics-inspired regularization term [69,70] constrains predicted temperature transitions to physically plausible rates consistent with land surface thermal inertia [71,72]. It is important to clarify the epistemological status of this constraint: the $L_{physics}$ term does not encode an explicit physical law such as a surface energy balance equation or a radiative transfer model. Rather, it encodes a physically motivated prior, that month-to-month LSTn transitions in continental arid environments are bounded by the thermal inertia of the land surface, and that transitions exceeding this bound are thermodynamically implausible regardless of data-driven fit. This use of domain knowledge as a soft regularization constraint is consistent with the broader “physics-guided machine learning” literature [73,74,75], in which the term encompasses any approach that embeds physically motivated constraints into the learning objective, whether derived from first principles, empirical physical bounds, or process understanding, without requiring the model to solve governing differential equations. Under this definition, the $L_{physics}$ term qualifies as physics-guided because its form and threshold are derived directly from the known thermal inertia properties of the land surface rather than from the training data alone.

$$L_{PG}=\alpha \cdot MSE+\beta \cdot L_{physics}$$

(17)

$$MSE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}(Y_i-{\widehat{Y}}_i{)}^2$$

(18)

$${\mathcal{L}}_{physics}={\displaystyle \frac{1}{n-1}}{\displaystyle \sum _{i=1}^{n-1}}\left[\right({\widehat{Y}}_{i+1}-{\widehat{Y}}_i)-(Y_{i+1}-Y_i){]}^2$$

(19)

To ensure reproducibility, Table 2 summarises the complete hyperparameter configuration used for all models in this study, where $\alpha =0.85$ and $\beta =0.15$ balance pointwise accuracy with temporal smoothness. This hybrid architecture directly addresses task mismatch, absence of long-range attention, and physically unconstrained training three common limitations in conventional LST forecasting studies.

3.4. Model Evaluation and Uncertainty Analysis

3.4.1. Performance Metrics

Model performance was quantified using three complementary metrics.

Coefficient of Determination (R²):

$$R^2=1-{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n}(Y_i-{\widehat{Y}}_i{)}^2}{{\displaystyle {\sum }_{i=1}^n}(Y_i-\bar{Y}{)}^2}}$$

(20)

where $Y_i$and ${\widehat{Y}}_i$are observed and predicted nighttime LST, respectively, and $\bar{Y}$is the mean of the observations. R² approaches 1 for accurate predictions, equals 0 when the model performs no better than the observed mean, and can be negative when predictions are worse than the mean.

Root Mean Square Error (RMSE):

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}}({\mathrm{Y}}_{\mathrm{i}}-{\widehat{\mathrm{Y}}}_{\mathrm{i}}{)}^2}$$

(21)

RMSE quantifies the magnitude of prediction errors in the original units (°C), providing a direct interpretation of forecast accuracy and penalizing larger errors more heavily.

Mean Absolute Error (MAE):

$$MAE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\mid Y_i-{\widehat{Y}}_i\mid$$

(22)

MAE provides a robust error measure that is less sensitive to outliers than RMSE. Seasonal performance was further evaluated across the four standard climate seasons - DJF, MAM, JJA, and SON - to assess model robustness across Afghanistan’s distinct seasonal regimes.

3.4.2. Forecast Uncertainty and Degradation Analysis

Month-by-month forecast errors were analysed to quantify error growth beyond the validated horizon. Linear regression of MAE values across validation years revealed a statistically significant degradation rate:

$$MA{E}_{degradation}=0.069°C/year(p=0.048)$$

(23)

This empirically established degradation rate underpins a three-tier confidence stratification framework for long-range projections:

High confidence (2025-2029): MAE = 1.34-1.61 °C

Medium confidence (2030-2034): MAE = 1.68-1.96 °C

Low confidence (2035-2040): MAE = 2.09-2.37 °C

3.4.3. Ablation Study and Component Contribution

To quantify the independent contribution of each PG-ST-LSTM component, four model variants were trained and evaluated under identical hyperparameter settings:

• Full PG ST LSTM (all components);
• Without physics-guided regularization (MSE only);
• Without Transformer encoder (BiLSTM retained);
• Without multi-driver inputs (LSTn and seasonal encoding only).

Component contributions were measured by the change in R² and RMSE upon removal. This ablation framework parallels the multi-condition uncertainty analysis employed in landslide susceptibility modeling, allowing robust quantification of which architectural elements drive performance improvements.

3.4.4. Trend Analysis and Trend-Constrained Projections

Long-term warming trends in nighttime LST were estimated using two complementary non-parametric methods. The Mann-Kendall test [76,77] assesses whether a monotonic trend exists without assuming a specific probability distribution, while Sen’s slope estimator [78] provides a robust, outlier-resistant estimate of trend magnitude:

$${\beta }_{\mathrm{Sen}}=\mathrm{median}\left({\displaystyle \frac{x_j-x_i}{j-i}}\right),\hspace{1em}\forall \hspace{0.17em}i<j$$

(24)

where $x_i$ and $x_j$ are LST observations at times $i$ and $j$. To account for serial autocorrelation in the LST time series, the modified Mann–Kendall test [79] was applied at the pixel level. Basin-scale Sen’s slopes were estimated at the pixel level and spatially averaged within each of Afghanistan’s five major watersheds. The area-weighted national mean slope ($\beta =+0.044$ °C/year) was applied to the 2001–2024 baseline to derive trend-constrained projections for 2040:

$${\mathrm{L}\mathrm{S}\mathrm{T}}_{n,2040}={\mathrm{L}\mathrm{S}\mathrm{T}}_{n,\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}}+\beta \cdot \mathsf{\Delta }t$$

(25)

where $\mathsf{\Delta }t$ represents the projection horizon from the end of the validation period (2024), expressed in monthly time steps, and $\beta$ denotes the estimated warming rate applied over the projection period. This trend-constrained approach avoids recursive error accumulation inherent in unconstrained autoregressive model rollout, providing a defensible framework for long-range climate assessment in data-scarce environments. We note that the approach assumes the historical warming rate persists through 2040; potential acceleration of warming due to climate feedbacks would render this a conservative lower-bound estimate.

4. Empirical Results

4.1. Preliminary Diagnostic Analysis

Before implementing the deep learning forecasting framework, a comprehensive suite of preliminary diagnostic analyses was conducted to ensure the statistical validity of the data, justify modeling assumptions, and establish the empirical foundation for subsequent model development. First, Table S1 (see the supplementary file) presents the descriptive statistics for the four study variables, including LSTn, NDVI, PREC, and SM, disaggregated by season and river basin across the study period. The seasonal panel reveals distinct thermal variability, with mean LSTn ranging from -8.481 °C in Winter to 13.564 °C in Summer, reflecting Afghanistan’s extreme continental climate characterized by severe cold seasons and intense summer heating. Summer exhibits the lowest PREC (mean=13.505 mm) and SM (mean=0.167), consistent with the region’s arid monsoonal pattern, while Spring records the highest mean PREC (35.883 mm), indicating snowmelt-driven moisture input. NDVI values remain consistently low across all seasons (0.026-0.143), confirming sparse vegetation cover typical of hyper-arid to semi-arid landscapes. At the basin level, the Helmand basin records the highest mean LSTn (16.296 °C) with the lowest NDVI variability (Std. Dev.=0.003), reflecting its predominantly desert landscape and thermal stability, whereas the Amu-Darya basin exhibits the coldest mean LSTn (-9.879 °C) owing to its high-altitude northern position. The Northern basin shows the widest LSTn range (-24.07 to 11.12 °C) and highest PREC variability (Std. Dev.=28.495 mm), suggesting strong orographic influences and inter-annual climatic instability.

Additionally, Table 3 presents the Pearson correlation coefficients between LSTn and the predictor variables across seasons (Panel A) and basins (Panel B), revealing spatially and seasonally heterogeneous relationships. The NDVI-LSTn correlation exhibits a notable seasonal reversal, shifting from a strong positive association in Winter ($r=0.820$) to a negative one in Summer ($r=-0.515$), reflecting vegetation’s dual thermal role, insulation during cold months and evapotranspirative cooling during the growing season. PREC and SM consistently show negative correlations with LSTn across most seasons, particularly in Spring ($r=-0.633$ and $r=-0.580$, respectively) and Autumn ($r=-0.494$ and $r=-0.451$), confirming their surface cooling effect through latent heat exchange and energy balance modulation. At the basin level, strong positive NDVI-LSTn correlations in the Amu-Darya ($r=0.811$), Harirude ($r=0.837$), and Northern ($r=0.899$) basins are consistent with the arid-region paradox, wherein vegetation co-occurs with higher soil moisture retention and greater nocturnal longwave re-emission, while the near-zero correlation in the Helmand basin ($r=-0.065$) reflects its hyper-arid character, where vegetation is too sparse to exert meaningful thermal influence.

Moreover, the normality of all study variables was assessed using three complementary tests, including Shapiro-Wilk (SW), Jarque-Bera (JB), and Kolmogorov-Smirnov (KS), applied across all five basins and four seasons (Figures S2–S5 in supplementary file). Results reveal that LSTn exhibits predominantly mixed or normal distributions across most basin-season combinations, particularly in Winter and Summer, with confirmed normality in the Amu-Darya Winter (SW: $p=0.3796$, JB: $p=0.9038$, KS: $p=0.9103$), Harirud Winter ($p=0.8483$, $p=0.6535$, $p=0.9514$), Helmand Winter ($p=0.2979$, $p=0.3823$, $p=0.8628$), and Kabul Winter ($p=0.1350$, $p=0.4680$, $p=0.4057$), suggesting that cold-season LSTn follows a Gaussian distribution. In contrast, PREC consistently deviates from normality across all basins and seasons ($p<0.05$ in all three tests), reflecting the highly skewed and episodic nature of precipitation in Afghanistan’s arid climate. NDVI and SM exhibit predominantly non-normal or mixed distributions, particularly in Spring and Autumn, owing to their pronounced seasonal dynamics and spatial heterogeneity. The pervasive non-normality of predictor variables across basins and seasons confirms that parametric assumptions underlying the LR technique are frequently violated, providing a compelling methodological justification for adopting distribution-free deep learning architectures, including TCN, LSTM, BiLSTM, and the targeted PG-ST-BiLSTM, which impose no distributional constraints on input features. Furthermore, the stationarity test results in Figure S6 reveal mixed outcomes across the ADF and KPSS tests, while the Phillips-Perron (PP) test consistently confirms stationarity across all variables, basins, and seasons, suggesting that apparent non-stationarity detected by ADF and KPSS likely reflects sensitivity to structural breaks and heteroscedasticity rather than true unit-root processes. Importantly, the confirmed stationarity under the more robust PP test, combined with the distribution-free and non-parametric nature of DL architectures, which do not require strict stationarity as a prerequisite, affirms that the dataset is well-suited for TCN, LSTM, BiLSTM, and the targeted PG-ST-BiLSTM models, which are inherently capable of adaptively learning both stationary and near-stationary temporal dynamics through their internal gating mechanisms, attention layers, and physics-guided regularization without imposing restrictive assumptions on the underlying data-generating process.

Finally, Table 4 reveals that all predictor variables, including PREC, SM, and NDVI, exhibit VIF values substantially below the commonly accepted threshold of 5, with a mean VIF of 1.669, confirming the absence of problematic multicollinearity among the selected predictors. These results validate that PREC, SM, and NDVI contribute independent, non-redundant predictive information to the modeling framework, justifying their simultaneous inclusion as input features in all forecasting models without inflating coefficient estimates or compromising model stability.

These preliminary diagnostic analyses justify the adoption of distribution-free DL architectures, TCN, LSTM, BiLSTM, and the targeted PG-ST-BiLSTM, which are uniquely suited to model non-linear feature interactions across multiple temporal scales while incorporating physics-guided regularization to enhance predictive accuracy and physical interpretability.

4.2. Baseline Analysis

4.1.1. Performance Analysis of DL Models

To achieve accurate and reliable LSTn forecasting, this study employs five predictive models, including LR as a statistical baseline alongside four deep learning architectures, namely TCN, LSTM, BiLSTM, and the proposed PG-ST-BiLSTM. Among these, PG-ST-BiLSTM constitutes the primary targeted model, trained under two distinct physics-guided loss configurations to assess both predictive performance and parameter sensitivity: the original configuration (α=0.85, β=0.15) and a robustness check configuration (α=0.90, β=0.10), where α governs the weight assigned to prediction accuracy through the Mean Squared Error (MSE) term and β controls the temporal smoothness regularisation penalty that enforces physically plausible transitions between consecutive predictions. All trained models were generated both in-sample historical fits (2001-2024) and out-of-sample future forecasts (2025-2040) via an autoregressive rolling prediction strategy, with evaluation metrics, coefficient of determination (R²), Root Mean Square Error (RMSE), and Mean Absolute Error (MAE), computed on the held-out test set and stored alongside all predictions for subsequent comparative analysis.

As shown in Table 5, the proposed PG-ST-BiLSTM achieves the highest overall performance, with the highest R² (0.977) and lowest MAE (1.534 °C). On RMSE, the BiLSTM baseline records a marginally lower value (2.017 °C vs 2.023 °C for PG-ST-BiLSTM), a difference of 0.006 °C that falls within the measurement uncertainty of the validation metrics and is therefore not meaningful in practical terms. Across the two metrics where differences are interpretively significant-R² and MAE the PG-ST-BiLSTM leads by clear margins: ΔR² = +0.011 and ΔMAE = -0.006 °C, over BiLSTM. The physics-guided framework’s superior R² confirms that it better explains the variance structure of nighttime LST dynamics, while the lower MAE confirms more accurate average predictions. The marginal RMSE difference reflects the known sensitivity of RMSE to a small number of large errors rather than systematic predictive superiority of BiLSTM.

4.1.2. LSTn Prediction using DL Models

In this section, LSTn is forecasted across five river basins in Afghanistan using a suite of DL architectures, TCN, LSTM, BiLSTM, and the proposed PG-ST-BiLSTM, trained on a multi-variate feature set comprising NDVI, PREC, and SM as predictor variables. Each model was trained on the historical period (2001-2024) under a seasonal stratification framework and subsequently used to generate out-of-sample forecasts for 2025-2040 using an autoregressive rolling prediction strategy.

For the 2025-2040 projection period, future values of the three environmental drivers (NDVI, PREC, SM) are not independently forecast. Instead, the autoregressive rolling strategy recycles the model’s own LSTn predictions back into the input window at each step, while NDVI, PREC, and SM are held at their 2001-2024 climatological monthly means for each basin - computed as the long-term average of each variable for each calendar month across the full observational record. This climatological mean driver assumption is consistent with the trend-stationary forcing scenario adopted in this study, wherein historical warming rates are assumed to persist through 2040 as described in the methodology and results sections, and is analogous to the fixed-driver approach commonly adopted in satellite-based long-range thermal projection studies operating in data-scarce environments where future driver scenarios are unavailable [80]. The sensitivity of projected LSTn trajectories to this assumption is acknowledged as a primary limitation in Section 5, where potential non-stationarity in driver-LSTn coupling strengths under sustained warming is discussed.

LR is additionally included as a statistical benchmark to contextualize the performance gains attributable to DL architectures (see Figure S7 in the supplementary file). Among the evaluated models, PG-ST-BiLSTM constitutes the primary proposed framework, distinguished by its integration of a Transformer encoder block, stacked Bidirectional LSTM layers, and a physics-guided loss function that jointly optimizes prediction accuracy and temporal smoothness, thereby embedding domain-specific physical constraints directly into the learning objective.

Figure 5 presents the historical fit (2001-2024) and future projections (2025-2040) of the targeted PG-ST-BiLSTM model under the primary physics-guided loss configuration (α=0.85, β=0.15), wherein 85% of the loss is attributed to prediction accuracy (MSE) and 15% to temporal smoothness regularisation. In the upper panel, the model demonstrates a strong capacity to reproduce the observed seasonal LSTn cycles across all five river basins, with the predicted series (solid lines) closely tracking the observed values (dotted lines) from approximately 2001 onwards, following an initial warm-up period during 2001–2004 that is attributable to the 12-month look-back sequence initialisation. The model faithfully preserves the inter-basin thermal hierarchy throughout the entire historical period, the Helmand basin consistently records the highest LSTn values (peak summer: ≈25 °C) owing to its hyper-arid desert character and minimal vegetation cover, the Kabul basin exhibits intermediate thermal amplitudes reflecting its semi-arid montane climate, while the Amu-Darya basin records the coldest winter minima (≈-25 °C) consistent with its high-altitude northern continental position, confirming that the physics-guided architecture successfully learns basin-specific thermal signatures. In the lower panel, the PG-ST-BiLSTM projects stable, physically coherent, and well-defined seasonal cycles across all five basins through 2040, maintaining the characteristic sinusoidal oscillation pattern without exhibiting autoregressive drift, amplitude decay, or inter-basin rank reversal across the 16-year forecast horizon. The gradual stabilization of seasonal amplitude variability observed from approximately 2030 onwards, particularly notable in the Helmand and Kabul basins, suggests that the temporal smoothness penalty (β=0.15) effectively suppresses unrealistic high-frequency fluctuations in long-range autoregressive predictions, demonstrating that the physics-guided loss formulation not only improves in-sample fit but also enhances the physical plausibility and temporal coherence of out-of-sample forecasts, thereby validating the methodological rationale for embedding domain-specific constraints directly into the model’s training objective.

In parallel with the targeted PG-ST-BiLSTM DL technique, LSTn dynamics across all five river basins were additionally forecasted using three established DL architectures, including TCN, LSTM, and BiLSTM, serving as comparative baseline models to benchmark the predictive performance of the primary framework and to validate the accuracy, reliability, and generalizability of the projected LSTn trajectories. The inclusion of these baseline architectures enables a rigorous and systematic comparative evaluation, wherein performance differences across models can be attributed to specific architectural innovations, namely the Transformer encoder, stacked bidirectional processing, and physics-guided loss regularization rather than to data characteristics or preprocessing choices, thereby strengthening the methodological credibility and scientific robustness of the forecasting framework presented in this study.

Figures S7-S10 in the supplementary file present the historical fits (2001-2024) and future projections (2025-2040) generated by the TCN, LSTM, BiLSTM, LR models across all five Afghan river basins, collectively serving as comparative evidence to contextualize and validate the superior performance of the proposed PG-ST-BiLSTM. All three baseline architectures successfully reproduce the broad seasonal cyclicity of LSTn during the historical period, preserving the inter-basin thermal hierarchy with the Helmand basin recording the highest LSTn amplitudes and the Amu-Darya basin the lowest, demonstrating their general capacity to learn dominant temporal patterns from the training data. However, notable differences in fit quality are apparent; the TCN exhibits the smoothest and most consistent historical reconstruction with the closest alignment to observed values, followed by LSTM. In the forecast panels, all three models project stable seasonal cycles through 2040 without systematic upward warming trends, reflecting the inherent limitation of purely data-driven architectures that lack physical constraints to enforce thermodynamically plausible long-range behaviour, a critical deficiency directly addressed by the physics-guided loss function (α=0.85, β=0.15) embedded in PG-ST-BiLSTM, which simultaneously optimises predictive accuracy and temporal smoothness to produce physically coherent projections that the baseline models structurally cannot achieve.

Moreover, the Mann-Kendall Sen’s slope analysis applied to the PG-ST-BiLSTM-projected LSTn values (2025-2040) is presented in Table 6, revealing a continued positive warming trend across all five river basins and all four seasons over the forecast horizon. All Sen’s slope values are positive, confirming a universal nighttime warming tendency in the projected period, with a grand mean of +0.0160 °C/yr across all basin-season combinations, suggesting a marginal acceleration of nighttime surface warming under the trend-stationary climatological assumption adopted in this study, wherein historical warming rates are assumed to persist through 2040. At the seasonal level, Winter exhibits the highest projected warming rate (+0.1174 °C/yr), followed by Summer (+0.0984 °C/yr), Autumn (+0.0545 °C/yr), and Spring (+0.0504 °C/yr), indicating that cold-season and peak-summer nighttime temperatures are projected to remain the most thermally vulnerable seasons through 2040. At the basin level, the Helmand basin records the highest projected mean Sen’s slope (+0.0796 °C/yr), consistent with its hyper-arid desert character, minimal vegetation buffering, and high nocturnal thermal emissivity, while the Amu-Darya basin exhibits the lowest projected warming rate (+0.0460 °C/yr), attributable to its orographic cooling influence and comparatively higher soil moisture retention. These spatially and seasonally heterogeneous warming projections collectively affirm that Afghanistan’s nighttime thermal environment is on a persistent warming trajectory through 2040, with Winter and Summer emerging as the most thermally vulnerable seasons, findings that carry significant implications for water resource management, agricultural planning, and ecosystem resilience across the country’s major river basins.

Furthermore, Figure 6 presents a comparative analysis of observed versus fitted nighttime LSTn values generated by the proposed PG-ST-BiLSTM model across six progressively varying physics-guided loss configurations, wherein the weighting parameter α is systematically increased from 0.70 to 0.95 while the temporal smoothness regularisation penalty β is correspondingly decreased from 0.30 to 0.05. The results demonstrate that predictive accuracy improves monotonically as α increases, reflecting the growing dominance of the Mean Squared Error (MSE) prediction term over the smoothness regularisation component, with the model achieving its highest R², lowest RMSE, and lowest MAE at α=0.95 — though the marginal performance gains beyond α=0.85 suggest that the primary configuration (α=0.85, β=0.15) already achieves a near-optimal balance between predictive fidelity and physical plausibility. A detailed spatial breakdown of observed versus fitted LSTn values using the proposed PG-ST-BiLSTM model is provided in the supplementary material, where Figures S11–S15 illustrate the model’s predictive performance disaggregated across each of the five Afghan river basins individually, enabling basin-specific assessment of fit quality and spatial generalisation capacity.

4.3. Robustness Checks

To further validate the robustness, predictive stability, and parameter sensitivity of the proposed PG-ST-BiLSTM framework, a secondary loss configuration was employed wherein the physics-guided loss weighting was adjusted to $\alpha =0.90$ and $\beta =0.10$, increasing the relative contribution of the Mean Squared Error (MSE) prediction accuracy term while proportionally reducing the temporal smoothness regularisation penalty, and the model was retrained and re-evaluated across all five river basins and four seasons for both the historical fitting period (2001-2024) and the future forecast horizon (2025-2040). This robustness check is designed to assess whether the model’s predictive outputs are sensitive to moderate perturbations in the physics-guided loss parameterization, and to determine whether the balance between data-driven accuracy and physical constraint enforcement remains stable within a plausible range of loss configurations. If the projected LSTn trajectories under the alternative configuration (α=0.90, β=0.10) closely reproduce those generated by the primary configuration (α=0.85, β=0.15), in terms of seasonal amplitude, inter-basin thermal hierarchy, and future warming trajectory, this constitutes strong empirical evidence that the PG-ST-BiLSTM architecture is not overly sensitive to the specific loss weighting scheme and that its forecasting behaviour is governed primarily by the learned spatiotemporal dynamics of the training data rather than by the precise calibration of the regularisation parameter, thereby confirming the scientific credibility and methodological reliability of the targeted DL framework.

As illustrated in Figure 7, the historical fit (2001–2024) and future projections (2025-2040) of the PG-ST-BiLSTM model under the robustness check configuration ($\alpha =0.90$, $\beta =0.10$), wherein a greater weight is assigned to prediction accuracy, and a reduced penalty is applied to temporal smoothness regularisation. The results demonstrate that the alternative loss configuration produces LSTn trajectories that are virtually indistinguishable from those generated by the primary configuration ($\alpha =0.85$, $\beta =0.15$), both in terms of seasonal cyclicity, inter-basin thermal hierarchy, and future warming behaviour, confirming that the model’s predictive outputs are robust to moderate perturbations in the physics-guided loss weighting scheme. In the upper panel, the predicted series (solid lines) closely track the observed values (dotted lines) across all five basins from approximately 2005 onwards, preserving the established thermal hierarchy wherein the Helmand basin consistently records the highest LSTn amplitudes and the Amu-Darya basin the lowest, with no meaningful degradation in fit quality relative to the primary configuration. In the lower panel, the future projections (2025-2040) maintain stable, physically coherent seasonal cycles across all basins without autoregressive drift or amplitude collapse, and the close correspondence between both loss configurations throughout the entire forecast horizon collectively constitutes strong empirical evidence that the PG-ST-BiLSTM architecture is parameter-stable, and that its forecasting behaviour is governed primarily by the learned spatiotemporal dynamics of the training data rather than by the precise calibration of the regularization weighting, thereby validating the scientific credibility and methodological reliability of the selected DL analytical framework in this study.

4.4. Discussion

This study set out to develop a robust deep learning framework for forecasting LSTn across Afghanistan’s heterogeneous river basins. The findings demonstrate that the proposed PG-ST-BiLSTM model significantly outperforms baseline statistical and DL methods, particularly by integrating a temporal smoothness constraint. A core contribution of this study is the explicit empirical justification for adopting deep learning over classical parametric models. The preliminary diagnostic analysis revealed that while LSTn can be normally distributed in cold seasons, consistent with the stable, Gaussian behaviour noted by Rani et al. [81] for nighttime LST predictor variables like PREC, they exhibit pervasive non-normality across all basins and seasons. As Harod et al. [82] also found when downscaling LST across India, the relationships between LST and its drivers are rarely linear or scale-invariant. Our results, showing that the Linear Regression baseline produced the weakest performance (R²=0.926), directly confirm that violating parametric assumptions degrades predictive fidelity. This empirically validates the necessity of distribution-free DL architectures (TCN, LSTM, BiLSTM, and PG-ST-BiLSTM) used in this study.

Additionally, the success of our PG-ST-BiLSTM’s physics-guided loss function (α=0.85, β=0.15) directly addresses a critical gap identified by Kim et al. [60], who called for “advancing modeling techniques” and integrating physical constraints into nighttime LST studies. The baseline TCN, LSTM, and BiLSTM models, while capturing seasonal cycles, produced purely data-driven projections and lacked a mechanism to enforce thermodynamically plausible long-range behaviour. By explicitly penalizing non-physical temporal discontinuities, our PG-ST-BiLSTM generated projections that maintained stable seasonal amplitudes without autoregressive drift through 2040 (Figure 5). The robustness check (α=0.90, β=0.10) further confirmed that this behaviour is not an artefact of hyperparameter tuning but a stable property of the architecture, providing a methodological template for embedding domain knowledge into DL as advocated by Gong et al. [83] for air temperature estimation.

Moreover, our empirical results reveal a complex, seasonally reversing relationship between LSTn and its drivers, which has direct implications for understanding arid-zone thermal dynamics. The seasonal reversal of the NDVI-LSTn correlation, from strongly positive in winter (r = 0.820) to negative in summer (r = -0.515), is a key finding. This supports the “arid-region paradox,” in which vegetation acts as an insulator during cold nights (trapping outgoing longwave radiation) and as a coolant via evapotranspiration during hot periods. This dual role is more nuanced than the general positive relation between built-up area and nighttime LST reported by Chakraborty et al. (2021) in an urban context, highlighting that in natural arid landscapes, vegetation’s thermal signal is seasonally contingent. The consistently negative correlations between PREC, SM, and LSTn across most seasons confirm their surface-cooling effect via latent heat exchange. This aligns with Muse et al. [84], who found that increased air moisture limits nighttime cooling in humid Florida. However, our findings invert that mechanism for arid Afghanistan: here, moisture input (PREC/SM) is the cooling agent, whereas in humid Florida, ambient air moisture limits cooling. This contrast underscores the importance of geographic and climatic context when interpreting LSTn drivers, a point also emphasized by Liang et al. [85] in their discussion of the decoupling of light and heat across different climate zones in China.

Furthermore, the baseline analysis provided a clear hierarchy of model performance based on R²: PG-ST-BiLSTM > LSTM = BiLSTM > TCN > LR. Among the baselines, TCN demonstrated competitive performance (R²=0.961) is consistent with its dilated causal convolutions, which provide a large receptive field. Notably, BiLSTM underperformed relative to unidirectional LSTM. This suggests that for seasonal LSTn forecasting, future context (the backward pass in a BiLSTM) may introduce noise rather than signal, possibly because the seasonal cycle is strongly dependent on a fixed historical look-back window, and bidirectional processing can lead to overfitting to interannual anomalies. Importantly, the proposed PG-ST-BiLSTM, despite being built upon a BiLSTM backbone, overcame this limitation by integrating a Transformer encoder and the physics-guided loss. The Transformer’s attention mechanism likely helped the model re-weight the importance of forward and backward temporal dependencies, while the smoothness penalty prevented overfitting to spurious patterns. This result echoes Harod et al. [82], who found that ResNet architectures performed well during the day but required separate nighttime models. Our unified PG-ST-BiLSTM, by contrast, successfully handled all seasons and basins, demonstrating the value of adaptive attention mechanisms over static architectural choices. The Mann-Kendall Sen’s slope analysis indicates persistent and widespread nighttime warming across Afghanistan (+0.0160 °C/yr, grand mean), with Winter (+0.1174 °C/yr) and Summer (+0.0984 °C/yr) emerging as the most vulnerable seasons. This finding is both methodologically and climatically significant. Methodologically, our projected warming trend in nighttime LST contrasts with the daytime cooling trends reported by Rani et al. [81]across most of India, using MODIS and AIRS data. While Rani et al. [81] found consistent nighttime warming (0.010-0.049 °C/yr) across most Indian climate zones, our rates are broadly comparable but slightly higher for winter. This difference may reflect the more extreme continental aridity and lack of moderating monsoonal influence in Afghanistan, leading to a stronger nocturnal radiative signal. Environmentally, the finding that winter and summer nights are warming fastest has profound implications. Cold-season warming can disrupt snowpack dynamics and water availability for spring snowmelt, a critical resource in the Helmand and Amu-Darya basins. Simultaneously, warmer summer nights reduce the diurnal temperature range, as reported by Muse et al. [84] and Logan et al.[86] What has been shown is a key driver of human thermal discomfort and heat stress, as the body has less opportunity to recover from daytime heat. The Helmand basin’s highest projected warming rate (+0.0796 °C/yr) is consistent with its hyper-arid, unvegetated character, which lacks the evapotranspirative buffering seen in other basins, confirming Logan et al.’s [87] finding that vegetation is the most important urban characteristic mitigating LST.

4.5. Limitations and Future Directions

Despite the robust performance of the proposed PG-ST-BiLSTM framework, three primary limitations should be acknowledged. First, the model does not explicitly account for anthropogenic drivers such as land use/land cover change, urbanization, or irrigation expansion, which may influence long-term LSTn trends beyond the climatological forcing considered here. Future research should integrate time-varying land surface modification layers to isolate natural versus anthropogenic warming signals. Second, the projected warming rates assume stationarity in the relationship between LSTn and its predictors (NDVI, PREC, SM) over the 2025-2040 horizon. However, feedback processes-such as drought-induced vegetation decline or soil moisture depletion-may alter these coupling strengths under sustained warming. Future studies should test for non-stationarity using rolling window correlation or regime-switching models. Third, uncertainty quantification is limited to parameter sensitivity (α, β weighting) and does not include propagation of input data errors, model structural uncertainty, or stochastic climate variability. Future work should adopt Bayesian deep learning or conformal prediction intervals to produce probabilistic LSTn forecasts with quantifiable confidence bounds.

4.6. Conclusions

This study introduced a novel Physics-Guided Transformer Bidirectional Long Short-Term Memory (PG-ST-BiLSTM) model for nighttime land surface temperature (LSTn) forecasting, using Afghanistan-a topographically complex and data-scarce region-as a challenging test bed. Through a rigorous suite of preliminary diagnostic analyses, baseline model comparisons, and robustness checks, the following key conclusions are drawn.

First, the preliminary analysis established a strong empirical and statistical foundation for deep learning-based LSTn forecasting. Descriptive statistics revealed distinct seasonal and basin-level thermal variability, with mean LSTn ranging from -8.481 °C in Winter to 13.564 °C in Summer, reflecting Afghanistan’s extreme continental climate. Pearson correlation analysis uncovered a seasonal reversal in the NDVI-LSTn relationship, shifting from positive in Winter to negative in Summer, confirming vegetation’s dual thermal role. Normality tests confirmed that precipitation (PREC) consistently deviates from normality across all basins and seasons, while stationarity tests under the robust Phillips-Perron (PP) test confirmed that all variables are stationary, collectively justifying the adoption of distribution-free deep learning architectures that impose no parametric constraints. Furthermore, variance inflation factor (VIF) analysis confirmed the absence of problematic multicollinearity among predictors (mean VIF = 1.669), validating the simultaneous inclusion of NDVI, PREC, and soil moisture (SM) as independent input features. Second, the baseline performance analysis demonstrated that the proposed PG-ST-BiLSTM model significantly outperforms both conventional deep learning architectures and statistical baselines. As shown in Table 4, PG-ST-BiLSTM achieved the highest predictive accuracy with an R² of 0.977, RMSE of 2.023 °C, and MAE of 1.534 °C, surpassing BiLSTM and LSTM (both R² = 0.966), TCN (R² = 0.961), and linear regression (R² = 0.926). The observed versus fitted analysis (Figure 6) confirmed strong alignment along the reference line (R² = 0.977), with a trend line of y = 0.98x + 0.21. These results conclusively demonstrate that the integration of a Transformer encoder for long-range attention, stacked BiLSTM layers for bidirectional sequence learning, and a physics-guided temporal smoothness regularizer provides a substantial and quantifiable improvement over purely data-driven architectures. Third, the projected LSTn trajectories (2025-2040) under the PG-ST-BiLSTM model reveal a persistent and spatially heterogeneous nighttime warming trend across Afghanistan. Mann-Kendall Sen’s slope analysis confirmed that all slope values are positive, with a grand mean of +0.0160 °C/year across all basin–season combinations. At the seasonal level, Winter exhibits the highest projected warming rate (+0.1174 °C/year), followed by Summer (+0.0984 °C/year), Autumn (+0.0545 °C/year), and Spring (+0.0504 °C/year), identifying Winter and Summer as the most thermally vulnerable seasons through 2040. At the basin level, the hyper-arid Helmand Basin records the highest projected warming rate (+0.0796 °C/year), attributable to its minimal vegetation buffering and high nocturnal thermal emissivity, while the Amu-Darya Basin exhibits the lowest rate (+0.0460 °C/year) due to orographic cooling and higher soil moisture retention. These spatially and seasonally heterogeneous projections affirm that Afghanistan’s nighttime thermal environment is on a persistent warming trajectory, with significant implications for water resource management, agricultural planning, and ecosystem resilience. Fourth, the robustness check using an alternative loss configuration (α = 0.90, β = 0.10) confirmed that the PG-ST-BiLSTM architecture is parameter-stable and not overly sensitive to the precise calibration of the physics-guided regularization weight. In summary, the PG-ST-BiLSTM model demonstrated generalizability, reproducibility, and physical coherence for nighttime LST forecasting in a data-scarce, mountainous region. The model successfully preserves inter-basin thermal hierarchies, maintains stable seasonal cycles without autoregressive drift over a 16-year forecast horizon, and produces physically plausible warming projections that are robust to parameter perturbations. These findings have direct and immediate relevance to climate adaptation, water resource planning, agricultural risk assessment, and cryosphere monitoring not only in Afghanistan but also in other data-scarce, topographically complex regions worldwide. Future work should extend this framework to integrate additional physical constraints (e.g., surface energy balance equations), incorporate ensemble climate forcing scenarios, and validate the model’s transferability to other arid and semi-arid mountain regions.

A potential concern with applying a ~79,853-parameter model to 156 training sequences per basin is overfitting. Several lines of evidence in this study collectively mitigate this concern. First, early stopping with a patience of 15 epochs was applied during training, halting optimization as soon as validation loss ceased to improve and thereby preventing the model from memorizing training-set patterns. Second, Dropout (rate = 0.3) was applied after each BiLSTM layer, providing stochastic regularization that reduces co-adaptation of learned features. Third, the physics-guided loss function (β = 0.15) acts as an additional regularize by penalizing high-frequency temporal fluctuations, constraining the solution space beyond what MSE alone would achieve. Fourth, the robustness check under the alternative loss configuration (α = 0.90, β = 0.10) produced virtually identical validation metrics and projection trajectories, confirming that the model’s behaviour is governed by learned data structure rather than precise hyperparameter calibration, a hallmark of well-generalizing rather than overfitted models. Fifth, the validation R² of 0.977 on the held-out 2015-2024 period a ten-year window the model never saw during training, is consistent with the training-period fit, with no evidence of the performance degradation that would be expected under significant overfitting. Taken together, these five lines of evidence provide reasonable assurance that the model generalizes beyond the training set within the context of this study. We nonetheless acknowledge that formal k-fold cross-validation was not conducted due to the limited sequence count, and that future work with larger datasets should adopt cross-validation as a standard practice.

4.7. Implications

The methodological framework demonstrated here is applicable to nighttime LST assessment in data-scarce mountain regions where satellite-derived time series constitute the primary long-term thermal record. However, direct transferability requires site-specific validation, as topographic complexity, cloud cover regimes, and available historical record lengths vary across regions.

Within Afghanistan, the projected warming carries distinct implications for water and agricultural security, though these remain conditional on the assumption of stationary predictor–response relationships. In the Helmand Basin (+1.27 °C by 2040), sustained nighttime warming may intensify evapotranspiration and agricultural water stress, particularly during the summer growing season. In the Amu Darya Basin (+0.736 °C by 2040), the projected warming approaches the freezing threshold, with potential but uncertain implications for Hindu Kush glacier-fed streamflow and downstream water availability. In the Kabul Basin (+1.24 °C by 2040), snowmelt-dependent urban water supplies may face altered timing and reduced volume of spring runoff, though quantitative hydrological modeling is required to confirm these effects.

Beyond Afghanistan, the PG-ST-BiLSTM’s demonstrated reliability under extreme topographic and data-scarcity conditions suggests potential transferability to analogous environments across Central Asia, the Hindu Kush-Himalaya, and the Andes. However, direct empirical validation in each target region remains necessary before operational deployment. Future applications should prioritize (i) retraining the model with locally available satellite products, (ii) calibrating physics-guided loss weights to regional thermal regimes, and (iii) validating against in-situ observations where accessible. Where conventional monitoring infrastructure is sparse or absent, physics-guided deep learning offers a reproducible and scalable alternative for long-term LSTn assessment, provided that continuous satellite records of sufficient length are available.