Forecasting Monthly Pan Evaporation in Al-Najaf, Iraq Using SARIMA and Month-Wise ARMA Models for Water-Resource Risk Assessment

2.1 Study area and data collection

This study is conducted in Al-Najaf, located in the central-southern region of Iraq (Figure 1), which is characterized by a semi-arid climate. This region experiences extreme temperature variations and high evaporation rates, particularly during the summer months, which directly impacts the local hydrological balance [6]. This study is based on 57 years of historical monthly pan evaporation records, from 1968 to 2024, with a total of 684 observations. These long-term data are essential for understanding climatic trends and cyclical fluctuations, which are necessary for constructing accurate stochastic models [18]. The data were obtained from the Iraqi Meteorological Organization and Seismology (IMOS), Ministry of Transport, measured using a standard Class A evaporation pan located at Al-Najaf meteorological station (station ID: 40670, latitude: 31.9500°N, longitude: 44.3166°E, elevation: 50 m a.s.l.) and expressed in mm month^-1.

The total monthly record of 684 observations was divided into two sets: a training set of 600 months (January 1968 to December 2017) for model estimation, and an independent validation set of 84 months (January 2018 to December 2024) to evaluate predictive performance. A total of thirty-four monthly observations (5.6% of the 600-month training record) were missing, all occurring strictly within the training period as short, isolated gaps not exceeding two consecutive months due to routine gauge maintenance and record archiving transitions. An inter-annual seasonal interpolation scheme—interpolating missing values across corresponding months of adjacent years—was applied for data imputation [26, 27]. All imputation statistics were derived exclusively from the training period. The 2018–2024 validation period contained complete observations with zero missing values.

To ensure no systematic bias was introduced into the imputation method, statistical properties before and after imputation were verified; the series mean shifted negligibly from 284.65 mm to 286.02 mm (<0.5%), and the standard deviation shifted from 177.71 mm to 178.49 mm (<0.45%), with formal two-sample difference and variance testing confirming distributional preservation without statistical bias (t = 0.131, p = 0.895, variance ratio F = 1.009, p = 0.913). Furthermore, a formal sensitivity analysis was performed comparing the inter-annual scheme with the Kalman seasonal state-space smoothing method. The inter-annual scheme demonstrated superior accuracy, achieving lower error metrics (RMSE, MAE, symmetric mean absolute percentage error (sMAPE) and MAPE), while overall model rankings, dynamical structures, and validation inferences (SARIMA vs. month-wise ARMA) remained identical in both approaches. Prior to modeling, the record was screened for abrupt structural breaks and temporal heterogeneity. A homogeneity test across the 1968–2024 record (using Pettitt’s test) showed a statistically significant shift (p < 0.05) centered around 2008. Official station data confirmed that the Najaf station maintained continuous operation at fixed geographic coordinates and elevation using a standard Class A pan (120.7 cm diameter, 25.4 cm depth) under the uniform observational guidelines of the World Meteorological Organization (WMO) [28], confirming that this shift represents a genuine climatic signal resulting from regional drought trends, and not an observational anomaly. This statistical shift around 2008 is corroborated by documented hydro-climatic records in central Iraq and Al-Najaf station, where the severe regional drought of 2007–2009 triggered sharp precipitation deficits, marked positive temperature anomalies, and extreme negative drought indices (RDI and SPI) across the Middle Euphrates [29-31].

To account for this non-stationarity, first-order differencing and seasonal differencing (d = 1, D = 1) were applied within the SARIMA modeling framework, effectively turning the series into a stationary process before parameter estimation.

Figure 1. Geographic location of the study area (Al-Najaf, Iraq)

2.2 Data transformation and stationarity

Time series stationarity is a fundamental requirement for statistical models [32]. To obtain this, the original series of monthly pan evaporation depth was examined utilizing the Augmented Dickey-Fuller (ADF) test [33]. The ADF test was used to assess the presence of a unit root in the conditional mean process. Variance heterogeneity and distributional skewness were assessed separately using descriptive diagnostics and motivated the Box-Cox transformation.

The optimal Box‑Cox transformation parameter was estimated by maximizing the log-likelihood function over a range of candidate (λ) values. To stabilize the variance, the standard piecewise Box-Cox transformation was applied (Eq. (1)) [34].

$Y_t=\left\{\begin{array}{cc}\frac{X_t^\lambda-1}{\lambda} & \text { if } \lambda \neq 0 \\ \ln \left(X_t\right) & \text { if } \lambda=0\end{array}\right\}$            (1)

where, Y_t is the transformed value at time t, X_t is the observed monthly pan evaporation depth at time t and λ is the Box-Cox transformation parameter.

For the SARIMA (600 months) training series, likelihood maximization yielded λ = 0.260. The transformed series was then assessed for stationarity [35], followed by Seasonal and, if justified, regular differencing.

For the month-wise ARMA, transformation requirements were assessed separately using a likelihood-based Box-Cox estimation and distributional diagnostics. Only December required a material transformation, with λ ≈ 0.01, close to the logarithmic transformation; the other monthly subseries were retained on their original scale because the transformation did not materially improve variance stability or the residual diagnostics. Then the twelve series were standardized independently by Eq. (2) [10].

$Z_{y, m}=\frac{X_{y, m}^*-\mu_m}{\sigma_m}$             (2)

where, Z_y,m: monthly standardized series, X^*_y,m: the transformed (or untransformed, if no transformation was applied) monthly value, μ_m and σ_m: monthly mean and standard deviation estimated from the training period only.

To strictly prevent data leakage and guarantee genuine out-of-sample evaluation, all preprocessing parameters were derived exclusively from the 600-month training record (1968–2017). Specifically, the Box-Cox parameter (λ = 0.260 for SARIMA; λ = 0.01 for December ARMA) and the month-specific standardization moments (µ_m and σ_m) were estimated solely on training data. During the 2018–2024 validation period, out-of-sample observations and their generated forecasts were back-transformed and unstandardized using these fixed training parameters without any re-estimation.

2.3 Mathematical formulation

SARIMA model is known as SARIMA (p,d,q) (P,D,Q)_S, whereas (p,d,q) refers to the non-seasonal part of the original series, (P,D,Q) refers to the seasonal part of the series and s is the season; this can be written as in (Eq. (3)).

$\varphi_p(B) \Phi_P\left(B^S\right)(1-B)^d\left(1-B^s\right)^D Y_t=\theta_q(B) \Theta_Q\left(B^S\right) \varepsilon_t$              (3)

where, p: the order of nonseasonal autoregressive (AR) model, q: the order of nonseasonal Moving average (MA) model, d: the number of regular differencing, P: the order of seasonal AR, D: number of seasonal differencing, Q: the order of seasonal MA, S: the length of season cycle, B: demonstrates the backward operators, Φ: seasonal AR coefficient, φ: non-seasonal AR coefficient, (1-B^s): seasonal differencing, (1-B)^d: non-seasonal differencing, θ: non-seasonal MA coefficient, and Θ: seasonal MA coefficient [36].

For comparison with SARIMA, the monthly record was decomposed into 12 calendar-month subseries, and an independent ARMA (p,q) model was fitted to each subseries across years. This approach captures month-specific interannual persistence but does not model dependence between adjacent calendar months; the mathematical model for each calendar month m is expressed in (Eq. (4)) [18]:

$Z_{m, y}=\sum_{j=1}^{p m} \emptyset_{m, j} * Z_{m, y-j}+\epsilon_{m, y}-\sum_{j=1}^{q m} \theta_{m, j} \epsilon_{m, y-j}$             (4)

where, Z_m,y is the standardized evaporation value for calendar month (m) in year (y); p_m and q_m denote the month-specific non-seasonal AR and MA orders, respectively; ϕ_m,j and θ_m,j represent the month-specific AR and MA parameters; and ɛ_m,y  is the independent white noise error term for month (m) in year (y).

The month-wise ARMA formulation is used here as a pragmatic seasonal-subseries benchmark rather than as a full periodic ARMA process. It captures interannual persistence within each calendar month but omits cross-month dependence; this limitation is explicitly considered when interpreting the comparison with SARIMA.

2.4 Model selection criteria

Information criteria were utilized to ensure parsimony and avoid overfitting. The Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) were utilized as shown in Eqs. (5) and (6):

$\mathrm{AIC}=-2 \ln (L)+2 k$              (5)

$\mathrm{BIC}=-2 \ln (L)+k \ln (n)$                (6)

where, L: the max likelihood function of the model, k: the number of estimated parameters in the model, and n: the number of observations in the model. The model that obtains the min value for both the AIC and BIC is considered the most suitable for representing the data [37].

2.5 Model identification

This is the first stage in the Box-Jenkins methodology, an overview of which is provided in Figure 2. Autocorrelation Functions (ACF) and Partial Autocorrelation Functions (PACF) are used to identify candidate time series models and their parameter orders (p,d,q) [17]. ACF and PACF plots were inspected using approximate 95% reference bounds of ±1.96/√N, where N represents the number of observations [38].

Figure 2. Box-Jenkins methodology plot

Table 1. Autocorrelation Functions (ACF) and Partial Autocorrelation Functions (PACF) criteria for identifying the order of Autoregressive Moving Average (ARMA) model

Candidate orders were selected from the observed cutoff/decay patterns and then compared by Information Criteria (AIC/BIC) to establish model optimum. Based on the statistical behavior of these functions, the methodological rules summarized in Table 1 were applied [19].

For the SARIMA framework, non-seasonal orders were systematically identified from correlation cutoff and decay patterns at early lags, while seasonal orders were inspected at seasonal cycle multiples (S = 12). To ensure a transparent and comprehensive model space, the search protocol evaluated candidate non-seasonal orders up to p,q ≤ 4 and seasonal orders P,Q ≤ 1 (with differencing fixed at d = 1, D = 1). For the month-wise ARMA framework, diagnostic identification was conducted independently on each calendar-month subseries across a parameter space bounded by p,q ≤ 2, establishing candidate AR and MA configurations based on the theoretical rules in Table 1 [18, 24].

2.6 Parameter estimation

This is the second stage in the Box-Jenkins methodology and aims to find the parameter values for the models identified in the first step. In this study, the model parameters were estimated utilizing MLE available in the SPSS software.

For the SARIMA model, the estimation algorithms were applied to the entire series. Stationarity and invertibility were verified by confirming that the roots of the AR and MA characteristic polynomials satisfied the required conditions (all inverse characteristic roots lie strictly within the unit circle) [17, 23]. In the month-wise ARMA model, however, the estimation process was performed for each month separately, with the parameters for each monthly series being extracted independently. This allows each month to express its own statistical characteristics without being bound by the assumption of parameter stability across the year [18, 24].

To confirm the significance of these parameters, t-tests were conducted for each parameter; parameters with a p-value < 0.05 were considered statistically significant. To formally verify the parsimony principle when eliminating statistically insignificant parameters, the Likelihood Ratio (LR) test was applied to evaluate constrained subset models against their unrestricted counterparts [39] (Eq. (7)):

$L R=-2\left(\ln L_{\text {restricted }}-\ln L_{\text {full }}\right) \sim \chi^2(d f)$               (7)

This parsimonious simplification ensured that there were no unnecessary coefficients that did not add explanatory value [11].

2.7 Diagnostic checking

After estimating the model parameters, the next stage is model diagnosis to assess the efficiency of data representation. This involves examining residuals, which are the differences between the actual values and the values predicted by the model. A model’s success depends on making the residuals behave like a “white noise” process. This means they must be independent, have a zero mean, and a constant variance [16, 40].

Residual independence was evaluated with the Ljung-Box Q test at lag h = 18 using degrees of freedom (df) adjusted for the number of fitted AR/MA parameters [41]. Residual normality was examined to assess the adequacy of Gaussian prediction-interval assumptions, using Q-Q plots and formal normality tests [42].

ACF and PACF were also plotted for the residuals, where the occurrence of all values within the confidence interval (±1.96/√N) indicates the absence of statistically significant correlations between them [19].

2.8 Performance evaluation metrics

To assess the efficiency of the nominated models in predicting monthly pan evaporation depth and to identify the model best able to represent the data, the statistical metrics demonstrated below were utilized to assess the performance of the SARIMA and month-wise ARMA models. These metrics are essential evaluation tools in hydrological research, as they provide indicators of the difference between the actual values (Y_t) and the values predicted by the models $\left(\widehat{Y}_{\mathrm{t}}\right)$ [23].

2.8.1 Root mean square error

Used to measure the overall error value, giving greater weight to large errors and penalizing extreme peak deviations, making it sensitive to outliers [43]. It is calculated utilizing Eq. (8):

$R M S E=\sqrt{\frac{1}{n} \sum_{i=1}^n\left(Y_i-\widehat{Y}_i\right)^2}$                (8)

2.8.2 Mean absolute error

Measures the average value of the actual errors on a continuous linear scale regardless of their direction, providing an unweighted baseline that avoids the quadratic inflation of RMSE [44]. It is calculated utilizing Eq. (9):

$M A E=\frac{1}{n} \sum_{i=1}^n\left|Y_i-\widehat{Y}_i\right|$                (9)

2.8.3 Mean absolute percentage error

Provides a relative standard that facilitates comparison between models [9]. It is calculated utilizing Eq. (10):

$M A P E=\left(\frac{1}{n} \sum_{i=1}^n\left|\frac{Y_i-\hat{Y}_i}{Y_i}\right|\right) \times 100$               (10)

2.8.4 Symmetric mean absolute percentage error

This measure assesses the percentage of relative error between observed and predicted values, giving equal weight to overestimates and underestimates [45]. It is calculated using Eq. (11):

$s M A P E=\frac{1}{n} \sum_{i=1}^n \frac{\left|Y_i-\hat{Y}_i\right|}{\left(\frac{\left|Y_i\right|+\left|\widehat{Y}_i\right|}{2}\right)} \times 100 $           (11)

2.8.5 Nash-Sutcliffe efficiency

This measure assesses the predictive power and ability of hydrological models to accommodate variance compared to the observed mean variability [46]. It is calculated utilizing Eq. (12).

$N S E=1-\frac{\sum_{i=1}^n\left(Y_i-\hat{Y}_i\right)^2}{\sum_{i=1}^n\left(Y_i-\bar{Y}\right)^2}$           (12)

2.8.6 Pearson’s coefficient of determination

It measures the strength of linear correlation and the consistency of structural patterns between observed and predicted values [47]. It is calculated using Eq. (13):

$Pearson$ $R^2=\left(\frac{\sum\left(Y_i-\bar{Y}\right)\left(\widehat{Y}_i-\mu_{\hat{Y}}\right)}{\sqrt{\sum\left(Y_i-\bar{Y}\right)^2 \sum\left(\widehat{Y}_i-\mu_{\hat{Y}}\right)^2}}\right)^2$              (13)

where, $\bar{Y}$ is the mean of observed values, $\widehat{Y}_i$ is the predicted value and $\mu_{\hat{Y}}$ is the mean of the predicted values. While Pearson $\mathrm{R}^2$ quantifies collinearity and structural alignment [47], it is evaluated alongside NSE to provide a strict hydrological benchmark that penalizes variance differences and systematic bias [46]. These metrics are applied to ensure comprehensive model accuracy. A model with the lower error values (RMSE, MAE, MAPE, sMAPE), higher Pearson $\mathrm{R}^2$, and NSE is considered optimal [48].

To evaluate whether the differences in forecast accuracy between the competing models are statistically significant, the Diebold-Mariano (DM) test was applied to the validation forecast error series under both squared-error and absolute-error loss functions [49], as shown in Eq. (14):

$D M=\frac{\bar{d}}{\sqrt{\frac{s_d^2}{N}}}$             (14)

where, $\bar{d}$: the mean loss differential ( $d_t=e_{1, t}^2-e_{2, t}^2$ for squared-error loss $d_t=\left|e_{1, t}\right|-\left|e_{2, t}\right|$ for absolute error loss), $d_t^2$: the sample variance of the loss differential series, $e_{1, t}$ and $e_{2, t}$ are the forecast errors of the competing models at time $t$, and $N$: the validation sample size.

All transformations, model orders, and parameters were strictly determined using the 1968–2017 training set only. Out-of-sample forecasts for the validation period (2018–2024) were generated using a fixed-origin multi-step forecasting protocol originating from December 2017. No future validation information was utilized at any stage of model selection, parameter estimation, or data preprocessing. Furthermore, forecasts generated on the Box-Cox transformed scale (λ = 0.260) were returned to the physical evaporation scale using the naive inverse Box-Cox transformation. Prediction limits were back-transformed consistently by applying the inverse transformation directly to the upper and lower bounds to preserve the specified coverage probability on the original physical scale prior to calculating evaluation metrics.

2.8.7 Seasonal naive benchmark model

To evaluate the practical forecasting skill and added value of the calibrated stochastic models, a seasonal naive approach was incorporated as a zero-parameter benchmark reference [50]. Consistent with the multi-step fixed-origin validation process originating from December 2017, the seasonal naive model projects future monthly pan evaporation by repeating the last observed annual seasonal cycle of the training series across the validation horizon $\left(\widehat{Y}_i=Y_{i-12}\right)$.

The research reveals that the evaporation time series for the Najaf area have a significant seasonal autocorrelation, making them naturally suitable for stochastic characterization utilizing the Box-Jenkins approach. During the 2018–2024 validation period, SARIMA and month-wise ARMA models produced comparable forecast errors, although the ranking depended on the metric: the month-specific model had a marginally lower RMSE (62.21 mm vs. 62.40 mm), while SARIMA achieved lower MAE, MAPE and sMAPE. Formal DM testing confirmed that this difference in forecast accuracy between the two models is statistically insignificant (p > 0.05).

A key feature of the month-wise ARMA model is that each calendar month is assigned its own interannual dependence structure, allowing heterogeneous month-specific dynamics. Its limitation is that it does not model dependence between adjacent months. SARIMA, in contrast, provides a single parsimonious seasonal structure and yielded lower MAE, MAPE and sMAPE in validation, motivating its use for the long-horizon reference forecast.

The long-horizon baseline forecast preserves pronounced seasonal peaks and troughs, while uncertainty widens with the forecast horizon. Crucially, these projections should be interpreted as a stochastic planning reference for atmospheric evaporation demand variability rather than a direct or deterministic assessment of future water losses, as actual volumetric losses remain contingent upon open-water surface dynamics and reservoir operational regimes. SARIMA and the month-specific model provide complementary statistical perspectives on temporal persistence and seasonal heterogeneity, offering regional water authorities a probabilistic baseline for water-resource risk assessment.

Despite the structured statistical workflow, the study has several limitations. First, it relies on observations from a single meteorological station in Najaf, limiting spatial generalizability. Second, the univariate models omit exogenous drivers such as temperature, wind speed, relative humidity, and radiation; these could be incorporated through SARIMAX or dynamic-regression models. Third, while the calibrated stochastic models offer structural transparency, they inherently assume linear temporal dependence, omitting complex nonlinear climatic interactions that data-driven machine-learning algorithms (e.g., LSTM and XGBoost) can capture. Future work should systematically benchmark these linear stochastic baselines against nonlinear machine-learning architectures under the identical time-aware protocol.