Forecasting Missouri riverflow with SARIMA models: A data-driven framework for adaptive water resourcemanagement

Introduction and Objectives

This study models and analyzes monthly flows of the Missouri River at Bismarck between 1954 and 2024 using a comprehensive dataset [15] . The water balance analysis explicitly incorporates precipitation, evaporation, upstream inflow, and tributary effects in the Upper Missouri River Basin [12] . The study applies ARIMA and SARIMA models to capture autocorrelation, long-term, and seasonal dynamics while assessing robustness through several statistical criteria and diagnostic tests [6,7,10,13] . The aim of overarching is to provide a reliable predictive framework to inform adaptive river management, conservation, and risk mitigation in a changing environment [11,16] .

Data and Methodology

Model fit is assessed via multiple metrics: Akaike Information Criterion (AIC), corrected AIC (AICc), Bayesian Information Criterion (BIC), and Nash-Sutcliffe Efficiency (NSE). Visual diagnostic plots with 95% credible intervals are used alongside quantitative metrics to provide a comprehensive performance evaluation, ensuring both statistical and physical validity [2,6] .

Description of Problem and Statistical Hypotheses

Suppose observed monthly water level at year 𝑡 follows 𝑦 𝑡 = 𝜇 𝑡 + 𝜖 𝑡 where 𝜇 𝑡 the average is acre-feet for monthly water level (month 𝑡), and 𝜖 𝑡 is random noise. Tests for stationarity include the Augmented Dickey-Fuller (ADF) and KPSS tests, applying hypotheses. The central statistical question concerns the stationarity and predictability of Missouri River water levels. We formally state: a) Null hypothesis (𝐇 𝐎 :): The series is stationary and thus suitable for conventional time series fore casting.

b) Alternative hypothesis (𝐇 𝛂 :): The series is nonstationary, requiring transformation or advanced models for accurate prediction [2,7] .

Model Structure and Testing Methodology

Stationarity is evaluated using the Augmented Dickey-Fuller (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests, with significance level p < 0.05 [4] . Detection of nonstationarity informs the use of SARIMA or Differencing approaches, accommodating trends and seasonality [5] . External drivers, especially climate change, are considered for their influence on underlying trends [13,14] .

Generation and Description of Data

Monthly flow data at USGS Site No. 06342500 (Bismarck) span January 1954-December 2024 (852 observations), with flow measured in cubic meters per second (m³/s) throughout [15] . Time series plots (not included: Figure 1) highlight typical ranges of 10,000-40,000 m³/s and a major flood event, spiking above 140,000 m³/s in 2010. Despite episodic surges, the long-term trend appears stable with no sustained increases or declines, suggesting underlying stationarity except during anomalies [1,12] . Figure 2 and Figure 3 illustrate the application of first-and second-order differencing to achieve stationarity. First, differencing reduces trends and stabilizes the mean (average); second differencing centers values around zero but may over-difference, potentially distorting real data structure.

Since the Augmented Dickey-Fuller (ADF) test confirms the original series is stationary, excessive differencing is avoided [15] . Figure 1 is the time series plot of the Missouri River's monthly flow from 1954 to 2024 highlights the river's dynamic hydrological behavior over seven decades. While the data generally fluctuates within a consistent range typically between 10,000 and 40,000 units there are periods of heightened volatility, especially during the 1970s and after 2010. A particularly striking feature is an extreme outlier around 2010, where flow values spike sharply above 140,000 units, likely corresponding to a major flood event. Despite these episodic surges, the long-term trend appears relatively stable, with no clear evidence of sustained increase or decline in flow levels. This pattern suggests that while the Missouri River is subject to extreme hydrological shocks, its overall flow regime remains stationary in the absence of such anomalies.

Stationarity Check and Differencing

Stationarity check in time series analysis involves determining whether the statistical properties of the data such as mean, variance, and auto covariance remain constant over time, which is crucial for reliable modeling and forecasting. If a series is non-stationary, it may have trends, seasonality, or changing variance, making it difficult to predict. Differencing is a common technique used to transform a non-stationary series into a stationary one by subtracting consecutive observations, which helps remove trends and stabilize the average (mean). The concept of a unit root relates to stationarity by testing whether the series has a root equal to one in its characteristic equation, which implies nonstationarity; if a unit root is present, the time series tends to have persistent shocks and evolving statistical properties. The Augmented Dickey-Fuller (ADF) test is used to determine whether a time series is stationary or contains a unit root, which would indicate non-stationarity, a crucial distinction in time series forecasting. The stationarity analysis of the Missouri River flow data confirms that the original time series is already stationary, making it appropriate for direct modeling with techniques like ARIMA. In this case, the ADF statistic for the original series was significantly negative (-5.4069) with a p-value below 0.01, leading to a rejection of the null hypothesis and confirming stationarity [2] . While differencing is commonly applied to achieve stationarity, additional differencing of this already stationary series led to even more extreme ADF statistics (e.g., -19.588), suggesting over-differencing. Overdifferencing can strip the data of meaningful structure and introduce unnecessary noise. Therefore, the strong ADF results indicate that the original flow series is statistic stable over time, with no need for differencing, and is ready for effective time series forecasting [7,14] . Figures 2 and3 illustrating the first-and second-differenced Missouri River flow data demonstrate the step-by-step transformation toward stationarity, a critical requirement for accurate time series modeling. The first-differenced series, created by subtracting each observation from the one before it, reduces long-term trends and partially stabilizes the mean, although some irregular fluctuations and spikes remain. The second-differenced graph, which applies differencing a second time, further flattens the series and centers the values more consistently on zero, indicating a higher degree of stationarity. This process shows how Differencing removes non-stationary elements like trends or seasonality, making the data suitable for models such as ARIMA that rely on stable statistical properties over time. However, because the original Missouri River flow series was already stationary based on the ADF test, the second differencing step may represent over-differencing, which could unnecessarily distort the underlying structure of the data.

Training and Testing Sets

Decompose the Training Data

In the analysis of Missouri River flow, decomposing the training data such as the monthly discharge from 1954 to 2010 plays a vital role in uncovering the underlying structure of the time series. This process separates the observed series into three key components: trend, which reflects long-term changes in river flow possibly due to climate variability or upstream management; seasonality, which captures predictable patterns like annual snowmelt or rainfall-driven fluctuations; and residuals, which represent irregular or random noise. Using decomposition tools such as decompose () or st l () in R allows analysts to visualize these components and assess the strength of seasonal and trend behavior in the Missouri River data. This insight informs model selection such as whether to include seasonal or differencing terms in an ARIMA model and helps detect structural breaks or outliers like major flood events. Ultimately, decomposition enhances understanding of river dynamics and improves the accuracy and interpretability of forecasting models used for water resource planning and risk management [10] . This decomposition is critical for understanding the distinct drivers of river flow and supports the development of more accurate and interpretable forecasting models for the Missouri River.

Outlier Detection and Adjustment

Outlier detection and adjustment are essential processes in time series analysis used to identify and mitigate the influence of anomalous data points that deviate significantly from the general pattern. These outliers, which may result from measurement errors, sudden events, or structural changes, can distort trend and seasonality estimates, leading to inaccurate models and misleading forecasts. Detection methods include statistical tests, residual analysis from fitted models, and decomposition techniques. Once identified, outliers can be adjusted through imputation, transformation, or model-based correction to preserve the integrity of the data while minimizing their disruptive impact. Proper handling of outliers enhances the reliability and accuracy of time series modeling and interpretation [1,8] .

Model specification and Results

In this study, model specification refers to the process of selecting and defining the appropriate time series model structure to capture the flow dynamics of the Missouri River at Bismarck. After confirming that the river's monthly discharge data from 1954 to 2024 was stationary using the Augmented Dickey-Fuller (ADF) test, a Seasonal Autoregressive Integrated Moving Average (SARIMA) model was specified to account for both short-term and seasonal patterns. The final model, SARIMA(1,0,1)(1,0,1) [12] , incorporates one non-seasonal autoregressive (AR) term, one non-seasonal moving average (MA) term, and seasonal components with a 12-month lag to reflect the river's annual hydrological cycle driven by snowmelt and precipitation. A Box-Cox transformation (λ ≈ 0.37) was applied to stabilize variance and meet normality assumptions [6,7,10] . This model structure was chosen based on statistical criteria such as AIC and BIC, as well as diagnostic checks to ensure a good fit.

The model specification effectively captures the Missouri River's seasonal flow variations and autocorrelation structure, enabling accurate forecasting that supports flood management, reservoir operations, and ecological planning [11,12] .

Box-Cox Transformation

The Box-Cox transformation is a statistical technique used to stabilize variance and make a time series more normally distributed, which is often a key assumption for many forecasting models. It transforms non-linear relationships into linear ones by applying a power transformation controlled by a parameter, lambda (λ). When λ = 1, the data remains unchanged; when λ = 0, the transformation becomes a natural logarithm. Intermediate values apply to various degrees of power transformations. This method is especially useful for time series with heteroscedasticity (non-constant variance) or skewed distributions, helping to improve the performance and accuracy of models like ARIMA or linear regression. The statement "Data shifted by 6558.035 for Box-Cox transformation" indicates that a constant value (6558.035) was added to the original time series data to ensure all values are positive before applying the Box-Cox transformation, which requires strictly positive inputs. This is a common preprocessing step when the data includes zero or negative values. The Box-Cox transformation was then applied to the shifted data, and the optimal lambda parameter (λ = 0.3716728) was estimated using the Guerrero method, which selects a λ that stabilizes the variance across different segments of the data [6] . This transformation helps in meeting the assumptions of normality and homoscedasticity for better model performance.

Fit Time Series Models on Box-Cox adjusted data

Fitting time series models on Box-Cox adjusted data involves applying forecasting models such as ARIMA, exponential smoothing, or others on the transformed dataset to improve model accuracy and validity. Since the Box-Cox transformation stabilizes variance and makes the data more normally distributed, it helps fulfill the assumptions underlying many time series models. By working with this adjusted data, models are less influenced by heteroscedasticity or skewed distributions, which can otherwise distort trend and seasonality estimates. After fitting the model, forecasts generated in the transformed space are usually inverted back to the original scale by applying the inverse Box-Cox transformation, ensuring interpretability of the results [7,10] .

ARIMA/SARIMA Methodology

The study uses ARIMA (p,d,q)(P,D,Q)[s] (Autoregressive Integrated Moving Average) notation. Here, p/d/q denotes order of the non-seasonal autoregressive, differencing, and moving average terms, respectively; P/D/Q are their seasonal counterparts, and s is season length (12 for monthly data) [7,10] . Model selection incorporates AIC, AICc, and BIC to balance complexity and fit [6] . Outliers are flagged with diagnostics and either transformed or dampened using robust modeling. The ARIMA (3, 1, 1) (0, 0, 2) [12] with drift model, applied to Box-Cox transformed Missouri River flow data (λ = 0.3716728), effectively captures both the short-term memory and seasonal characteristics of the river's discharge patterns.

The non-seasonal AR(3) and MA(1) components indicate that each month's river flow is influenced by the preceding three months of flow and one prior error term, reflecting the river's persistence and variability. The seasonal component with two seasonal moving average (SMA) terms and a 12-month cycle accounts for predictable annual hydrological behaviors, such as spring snowmelt or summer rainfall impacts. The inclusion of a small negative drift (-0.0519) suggests a gradual longterm decrease in river flow, potentially linked to environmental changes or upstream water management. The model's strong statistical fit indicated by a relatively low residual variance (σ² = 80.85), a high log-likelihood (-2466.99), and favorable AIC and BIC scores demonstrates its capacity to model the complexity of the Missouri River's behavior. This makes it a robust tool for forecasting river flow and supporting water resource planning across seasonal and inter annual time scales [7,10] .

Model Fit Statistics

Model selection by these criteria (particularly AIC/BIC) confirms good fit, while the relatively low residual variance indicates successful variance reduction. The table 3 shows the fit statistics from the ARIMA (3, 1, 1) (0, 0, 2) [12] with drift model offer valuable insight into how well the model captures the river's discharge dynamics. [8] .

Residual Diagnostics

Model adequacy is evaluated by residual plots, autocorrelation (ACF/PACF), the Ljung-Box test, and normality (Q-Q plot): Model diagnostics is the process of evaluating whether a fitted time series model, such as ARIMA, adequately captures the structure of the data and meets key statistical assumptions. It involves analyzing the residuals (the differences between observed and predicted values) to ensure they resemble white noise, that is, they should be random, normally distributed, have constant variance, and show no autocorrelation. Common diagnostic tools include residual plots, ACF/PACF plots of residuals, the Ljung-Box test for autocorrelation, and Q-Q plots for assessing normality. If these diagnostics indicate that assumptions are met, the model is considered reliable for forecasting and inference; otherwise, the model may need refinement [10] .

Ljung-Box Test Results (Residual Diagnostics)

While SARIMA addresses temporal structure more effectively than ETS (Exponential Smoothing) and TBATS (exponential smoothing with Box-Cox, ARMA terms, trend, and seasonality), significant autocorrelation remains, implying that additional lags, exogenous factors, or hybrid models could further improve performance. ETS and TBATS are less effective for the complex, bursty seasonality in this river data, leading to over-/under-estimation, as shown by their poor residuals. [12] 42.757 18 0.00087 Significant autocorrelation in residuals ARIMA(1,1,2)(2,0,1) [12] 36 (2, 0, 1) [12] Model Examines error patterns for signs of autocorrelation and normality. Some autocorrelation remains.

The residual diagnostic plots for the ARIMA(1,1,2)(2,0,1) [12] , ETS(A,N,A), and TBATS models highlight key differences in how well each model captures the data's structure. The ARIMA model shows relatively stable and normally distributed residuals, suggesting it captures much of the underlying trend and seasonality; however, minor spikes in its ACF plot indicate that some autocorrelation remains. The ETS(A,N,A) model demonstrates a tighter distribution of residuals in the time plot, but its ACF plot reveals more pronounced autocorrelation across several lags, signaling that it misses important temporal patterns in the river's flow.

Meanwhile, the TBATS model, despite being tailored for complex seasonal data, shows the most volatile residuals with a wider spread and significant autocorrelation, suggesting an overfit or structural mismatch. Overall, while none of the models fully satisfy the white noise assumption, ARIMA (1, 1, 2) (2, 0, 1) [12] offers the best residual diagnostics of the three, though additional model tuning or hybrid approaches may be needed to improve forecast accuracy for the Missouri River. Accurate test-period forecasting supports more reliable water resource planning, flood risk management, and environmental decision-making along the Missouri River [9.10] . Accurate evaluation on the test set ensures the chosen model is not only well-fitted to historical data but also effective for real-world decision-making [10,11] . through 2025, the model captures the full range of observed trends, seasonal cycles, and anomalies such as floods or droughts. This complete fit enhances the model's ability to produce accurate and robust future forecasts, which are essential for effective water resource management, infrastructure planning, and environmental decision-making related to the Missouri River [7] .

Table 6: Estimated Parameters for SARIMA (1, 0, 1) (1, 0, 1) [12] Model (Full Dataset) Tables 6 and 7shows the SARIMA(1,0,1)(1,0,1) [12] Based on selection criteria and diagnostics, the preferred SARIMA model for Missouri River monthly flow forecasting is: ARIMA (1, 0, 1) (1, 0, 1)12 with non-zero mean [12] The estimated parameters include: The residual variance isσ 2 = 68.31. This model captures both short-term autocorrelation and annual seasonal cycles, providing a balanced framework for hydro-meteorological time series forecasting [12] . The model equation using the backshift operator B is ARIMA (1, 0, 1) (1, 0, 1) [12] ∅ p (B)Φ P (B s )y t = θ q (B)Θ Q (B s )ε t y t = ∇ d ∇ s D ε t , ε t ~WN(0, σ 2 ) AR and MA Polynomials in the lag operator are respectively given by

The seasonal AR (SAR) and the seasonal MA (SMA) polynomials in the lag operator are respectively given by

(1, 0, 1)(1, 0, 1) [12] = (p, d, q) (P, D, Q)

Explicit solution for y t: y t -0.7044y t-1 +0.8315y t-12 +0.5856y t-13 =106.52+ε t +0.2880ε t-1 -0.6397ε t-12 [7] 4.7 Forecasted Monthly Missouri River Flow (ac-ft) with 95% Confidence Intervals

The forecasted monthly Missouri River flow from September 2024 to August 2025 shows a clear seasonal pattern, with flows gradually increasing from early fall to a peak in late spring (May 2025 at 23,175 ac-ft) before declining through the summer months. This trend aligns with typical hydrological cycles influenced by snowmelt and precipitation, with the highest flow occurring during spring. Confidence intervals are widest in high-flow months, particularly in March through May, indicating greater uncertainty and potential variability due to climatic factors. In contrast, fall and winter months exhibit more moderate flows and relatively narrower intervals. Overall, the projections highlight the need for proactive water resource planning, particularly in spring when both flood potential and forecast uncertainty are highest [12] . Final 13: Forecast: Missouri River Flow (2024-2025) Using SARIMA (1, 0, 1) (1, 0, 1) [12] Figure 13 displays a time series forecast of Missouri River flows using the best SARIMA model. The historical data shows high variability and several extreme flow events (e.g., around 1993 and 2011), while the forecast (in blue) for the near future indicates relatively stable flows with narrow confidence intervals, suggesting the model predicts moderate and consistent river behavior ahead.

Discussion

The Together, these parameters indicate that Missouri River flow exhibits both enduring seasonal patterns and autoregressive behavior, while also capturing transient shocks. These flow dynamics affect river ecosystem health by influencing water quality, aquatic habitat conditions, fish population dynamics, and shoreline erosion processes. Understanding these patterns through SARIMA aids in adaptive ecosystem and resource management [11,16] .

Conclusion

This study demonstrates that SARIMA models, particularly those enhanced with Box-Cox transformations, provide a robust and interpretable framework for forecasting monthly Missouri River flows. These models effectively capture both short-term autocorrelation and long-term seasonal dynamics, supporting adaptive decision-making in water resource management, reservoir operations, flood risk mitigation, and ecosystem planning. Among the models evaluated, the SARIMA (1, 0, 1) (1,0,1) [12] configuration offered the most balanced performance, combining low residual variance with strong fit statistics.

Although residual diagnostics, such as the Ljung-Box test, indicated some remaining autocorrelation, SARIMA models clearly outperformed alternatives like ETS and TBATS, particularly in capturing the complex seasonal patterns characteristic of river discharge data. Despite their strengths, SARIMA models depend heavily on historical patterns and may be less effective in forecasting rare or unprecedented hydrological events, such as extreme flooding. Therefore, continued model refinement is necessary. Future research should explore hybrid approaches that incorporate exogenous climate variables (e.g., temperature, snowpack, and precipitation forecasts) as well as non-linear or machine learning methods to better account for structural changes and reduce residual autocorrelation. From a policy and management standpoint, the use of advanced time series methods such as SARIMA strengthens the reliability of hydrological forecasting, providing a datadriven foundation for proactive and resilient water resource planning. Regular model updates and integration with climate projections will be critical to ensuring that river management strategies remain responsive to evolving environmental challenges [7,10] .

Acknowledgement

AcknowledgementsThis work is supported in part by the National Science Foundation grant no. HRD 1839895.

References

1. Abraham B, Chuang A. Outlier detection and time-series modelling. Technometrics. 1989 May;31(2):241-248.
2. Asare B, Sah M, Hona RK. Alternative methods of regular and singular perturbation problems. Appl Math. 2024;15:687-708. DOI: 10.4236/am.2024.1510040.
3. Asare B. Epsilon variations and their impact on solution accuracy: A matched asymptotic approach. Int J Stat Appl Math. 2024;9(6):34-44. DOI: 10.22271/maths.2024.v9.i6a.1883.
4. Asare B. Enhancing precision in nonlinear oscillator analysis: A comparative study of multiple-scale analysis and strained coordinates. Int J Stat Appl Math. 2025;10(1):128-141. Submitted. DOI: 0.22271/maths.2025.v10.i1b.1956.
5. Asare B, Rhonda M. Nonparametric tests for the umbrella alternative in a mixed design. Int J Stat Appl Math. 2025;10(6):64-70. DOI: 10.22271/maths.2025.v10.i6a.2060.
6. Akaike H. On the use of a linear model for the identification of feedback systems. Ann Inst Stat Math. 1968;20:425-439.
7. Box GEP, Jenkins GM, Reinsel GC. Time series analysis: Forecasting and control. 3rd ed. Englewood Cliffs, NJ: Prentice-Hall; 1994. 600 p.
8. Chatfield C, Prothero DL. Box-Jenkins seasonal forecasting: Problems in a case study (with discussion). J R Stat Soc Ser A. 1973;136:295-352.
9. Kumar P, Coulibaly P. River flow forecasting using artificial neural networks. J Hydrol Eng. 2009 MarApr;14(2):144-155.
10. Papacharalampous G, Tyralis H. Time series forecasting using SARIMA models: A review. Int J Forecast. 2018;34:610-632.
11. Salas JD, Obeysekera J, Vogel RM. Stochastic hydrology for the next millennium. J Hydrol Eng. 2010 Aug;15(8):581-587.
12. Schaefli B, Hingray B, Musy A. Short-term hydrological forecasting for flood-risk management. J Hydrol. 2005 Mar 10;309(1-4):328-342.
13. Spencer DE. Developing a Bayesian vector autoregression forecasting model. Int J Forecast. 1993;9:407-421.
14. Sutcliffe A. Time-series forecasting using fractional differencing. J Forecast. 1994;13:383-393.
15. Young PC. Data-based mechanistic modeling of environmental, ecological, economic and engineering systems. Environ Model Softw. 2011 Apr;26(4):331-345.