Article

September 30, 2026

Decision support in contaminated site remediation with machine learning algorithms

Decision Support in Contaminated Area Remediation with Machine Learning Algorithms

Deivid Lucas dos Santos Migueleti; Aimée Shirozono

DOI: 10.22167/2675-6528-202602835

Article derived from a Final Course Work (TCC), with content based on the student’s original work and adapted to the editorial format of the E&S Magazine with the support of the ResumeAI tool, an artificial intelligence solution developed by Instituto Pecege for textual synthesis and organization.

Summary

The remediation of areas contaminated by hydrocarbons in aquifers subject to seasonal water level fluctuations demands quantitative instruments for planning and optimization, whose efficiency is limited by the trapping of petroleum derivatives. The work demonstrated the applicability of statistical machine learning models in predicting time series of hydraulic head in an area contaminated by aviation kerosene in the Paulínia-SP petrochemical hub, as evidence for decision support in remediation. Monthly hydraulic head records from monitoring wells and monthly accumulated rainfall over a 10-year period were used. Time series behavior analysis, variable engineering, and correlation analysis aided in the estimation of SARIMA and SARIMAX models. Models with adequate predictive capacity were selected and adjusted based on parsimony criteria and residual diagnostics. Time series modeling with exogenous variables improved model performance compared to SARIMA models. The model with 3-month accumulated precipitation showed a mean absolute error (MAE) of 0.1041 meters for a 12-month forecast, which was considered satisfactory. This constituted a relevant proof of concept to support remediation planning and long-term management of hydrocarbon-contaminated areas.

Keywords: Groundwater; Contamination; LNAPL; Rainfall; Soil.

1. Introduction

Groundwater represents 99% of the planet’s freshwater reserves, being a fundamental strategic resource for the climate adaptation of societies. In this context, public policies aimed at the protection and recovery of underground springs are of extreme relevance (UNESCO, 2022).

Soil and groundwater contamination, however, can severely compromise the use of this vital resource. In the State of São Paulo, the Information System for Contaminated and Rehabilitated Areas (SIACR), maintained by the Environmental Company of the State of São Paulo (CETESB), registered, in October 2025, 7,253 contaminated areas, of which 4,846, approximately 66.8%, were affected by petroleum-derived automotive fuels (CETESB, 2025).

Petroleum derivatives present in the subsurface in their pure form are known as LNAPL, an acronym for Light Non Aqueous Phase Liquid, or immiscible non-aqueous phase liquid. These complex mixtures can contain toxic substances, such as benzene and naphthalene, which contaminate groundwater and pose risks to human health (Alden et al. 2024).

The rehabilitation process of a contaminated area involves intervention steps to mitigate or reduce the risks of exposure to chemical substances of interest (CETESB, 2017). Information on the distribution of contaminants and their interaction with the physical environment and receptors is consolidated in the Conceptual Environmental Model (CEM) (Alden et al. 2024).

Remediation measures use technologies to treat the affected environment, removing contaminants or restricting their mobility (containment) (CETESB, 2017). The São Paulo State Decree number 59.263/2013 prioritizes measures that promote the removal and reduction of contaminant mass (ALESP, 2013). The most common remediation techniques for areas contaminated by petroleum derivatives are based on the hydraulic removal of LNAPL from the aquifer, including multiphase extraction, pumping and treatment, free-phase recovery, and hydraulic barrier (ITRC, 2023; CETESB, 2017).

Despite its wide application, the efficiency of hydraulic removal techniques is limited by several factors and must be optimized as the area’s MCA is updated (ITRC, 2023). One of the phenomena that restricts LNAPL recovery is trapping, where LNAPL becomes imprisoned in soil pores below the water table due to differences in capillary forces and interfacial tension between oil and water (Isler et al. 2018; Teramoto et al., 2020). This phenomenon was observed by Pede (2009) in an area of the Paulínia petrochemical complex.

In regions with well-defined rainfall regimes, LNAPL trapping can be intensified by seasonal water table fluctuation in response to precipitation (Teramoto et al., 2020). Climate change and extreme events, such as prolonged droughts and torrential rains, add complexity to the management of contaminated areas (UNESCO, 2022). Rainfall well below average intensifies trapping, hindering remediation by hydraulic methods.

Although mathematical and computational models for simulating and predicting aquifer water levels have existed for a long time, their application is restricted due to complexity and the need for specialized resources. Previous research has used Artificial Neural Networks, Support Vector Machines, and statistical models such as ARMA and ARIMA to predict water levels, some with the ability to anticipate up to four months in advance (Tao et al., 2022; Nayak et al. 2006). More recently, models based on the Transformer architecture have demonstrated superior predictive capability (Sun et al. 2024).

The need for robust and more accessible machine learning models, capable of predicting water level, is evident. Such models, combined with knowledge about the vertical distribution of LNAPL, can provide a crucial line of evidence for decision-making and negotiation with environmental agencies regarding the discontinuation of remediation efforts in contaminated areas where LNAPL trapping has occurred. The objective of this work was to demonstrate the applicability of at least one machine learning model in predicting water level time series in an LNAPL-contaminated area as a subsidy for decision-making about the continuity of remediation through hydraulic removal techniques.

2. Material and Methods

This study was characterized as a quantitative research, with a methodological approach focused on time series analysis. The objective was to demonstrate the applicability of machine learning models in predicting the water level in an area contaminated by LNAPL, providing support for decision-making on the continuity of remediation by hydraulic techniques. The research was carried out in an area contaminated by aviation kerosene in the petrochemical complex of Paulínia, São Paulo, where the leak was stopped in 2006 and the LNAPL trapping phenomenon was previously observed by Pede (2009).

The unit of analysis consisted of monthly hydraulic load records from 118 monitoring wells (January 2013 to September 2025) and 180 monthly accumulated rainfall records (January 2011 to December 2025). The data were collected from monitoring wells and a surface weather station at the same industrial plant. The use of this company-owned data was authorized by a letter of agreement. The raw data were initially processed in Microsoft PowerBI Desktop for cleaning and treatment, and the hydraulic load was calculated by subtracting the water level from the altimetric elevation of the wellhead.

The subsequent steps of treatment, filling in missing information, and data preparation were conducted in the Python language, in the Google Colaboratory environment (Google, 2026), using libraries such as Pandas, NumPy, Matplotlib, Seaborn, Scipy, Scikit-learn, Statsmodels, and pmdarima. To fill gaps in the time series, the `interpolate` function from Pandas was employed. For model estimation, the period from January 2016 to September 2025 was used, excluding the water crisis period of 2014-2015. From the original variables, six lagged hydraulic load variables (one to six months) and five precipitation variables based on the Standardized Precipitation Index (SPI) for one, two, three, six, and twelve-month scales were calculated, using the SPEI package (Vonk, 2025), according to McKee et al. (1993). Accumulated precipitation variables were created using the `shift` function from Pandas.

An exploratory clustering step of the hydraulic load records was performed, using the DBScan algorithm from the Scikit-Learn library (Ester et al., 1996), after standardizing the values with StandardScaler. The `eps` and `min_samples` parameters of DBScan were adjusted based on a k-neighbors graph. Exploratory data analysis included descriptive analyses, box-plot graphs, line and bar graphs, and a Pearson correlation matrix, visualized through a heatmap generated with the Seaborn package.

The decomposition of the hydraulic load time series was performed using the `seasonal_decompose` function from the Statsmodels library (Hyndman and Athanasopoulos, 2021). Autocorrelation and partial autocorrelation (ACF and PACF) were evaluated with the `plot_acf` and `plot_pacf` functions from the same library. The stationarity of the time series was verified through the Augmented Dickey-Fuller (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests, using the `adfuller` and `kpss` functions from the `tsa` module of the Statsmodels library, at the five percent significance level (Morettin and Toloi, 2018). The hydraulic load series required two differencing steps to achieve stationarity.

For the estimation of the models, the `auto_arima` function from the pmdarima package was used, which minimizes the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) (Morettin and Toloi, 2018). Autoregressive Integrated Moving Average (ARIMA) models and their seasonal extensions (SARIMA), as well as SARIMAX models, which incorporate explanatory exogenous variables (Nelson, 1973), were estimated. The adjustment of the selected models considered only exogenous variables with statistically significant coefficients. Performance was evaluated by predictions of hydraulic load for twelve months, using the functions `get_forecast`, `predicted_mean`, and `conf_int`. The error metrics used were Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) (Hyndman and Athanasopoulos, 2021).

The diagnosis of the residuals was performed to verify their representativeness, including the generation of standardized residual plots, histograms, Q-Q plots, and correlograms. Formal verification of the residual characteristics included the Ljung-Box test for autocorrelation, and the Jarque-Bera, Shapiro-Wilk, and Heteroskedasticity tests for adherence to normality and homoscedasticity (Fávero and Belfiore, 2024).

3. Results and Discussion

The analysis of data from the area contaminated by aviation kerosene in Paulínia-SP revealed hydrological and pluviometric patterns that directly influence the behavior of the aquifer’s hydraulic head, providing essential subsidies for the application of predictive models. Initially, the exploration of hydraulic head records from monitoring wells allowed the identification of spatial groupings, although detailed cluster information was not incorporated into subsequent analyses to focus on the general behavior of the water table in the area. The spatial distribution of monitoring wells, as observed, showed a grouping that followed the direction of groundwater flow, from northeast to southwest, with some wells classified as “outliers” due to discrepant measurements.

The time series of hydraulic head, composed of 106 water level records after interpolating 9.40% of missing data, showed a median of 589.092 meters and a mean of 589.093 meters. Linear interpolation was applied to fill gaps, a recognized methodological limitation, but whose impact was considered limited by the nature of the aquifer’s response to precipitation and the short duration of measurement failures. The monthly distribution of hydraulic head revealed a clear seasonal pattern, with peaks in March and April and minima in October and November, indicating seasonal influence on the water level.

The monthly accumulated precipitation, with 117 records, showed a median of 75.000 millimeters and an average of 99.150 millimeters, with a slight right skew. The monthly precipitation distribution confirmed the expected seasonal pattern for the region, characterized by rainy summers and dry winters. This seasonality is a crucial factor for understanding the aquifer recharge dynamics and, consequently, the hydraulic head fluctuations, which are fundamental for planning the remediation of areas contaminated by LNAPL.

Correlation between variables

The evaluation of the correlation between the original and constructed variables, including hydraulic load, monthly rainfall, lagged hydraulic load variables, and standardized precipitation indices (SPI), revealed significant relationships. A 41% correlation was observed between the hydraulic load values and the 6-month SPI (spi6), 41% with the accumulated precipitation over 2 months (pluv_acum2), and 54% with the accumulated precipitation over 3 months (pluv_acum3). These results are consistent with the understanding that the aquifer’s response to precipitation is not immediate but rather delayed, due to the processes of infiltration and water recharge in the soil.

The analysis of hydraulic load time series together with monthly accumulated rainfall, the 6-month SPI, and 3-month accumulated rainfall, visually demonstrated that the hydraulic load tends to increase in periods of higher rainfall, confirming the dependency relationship. This observation is vital for modeling, as the inclusion of exogenous variables that capture this temporal lag can significantly enhance the predictive capacity of the models, aligning with the literature that highlights the importance of rainfall in aquifer dynamics (Teramoto et al., 2020).

Analysis of the decomposition of the response variable

The decomposition of the hydraulic load time series, performed by the additive method, clearly identified the presence of trend, seasonality, and residual components. The observed trend in hydraulic load over time reflects long-term changes in the hydrogeological system, while seasonality captures the annual patterns of water level oscillation, strongly influenced by the region’s rainfall regime. The identification of these components is a fundamental step for time series model estimation, as it allows models to capture and explain these behavioral patterns (Hyndman and Athanasopoulos, 2021).

The Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots of the response variable confirmed the existence of trend and seasonality in the data. The peaks in the PACF indicated that the time series of hydraulic load could be modeled as an autoregressive process of order 2 or 3. This information was crucial for choosing the initial parameters of the `auto_arima` function from the `pmdarima` package, guiding the structure of the SARIMA and SARIMAX models to be estimated.

Time series stationarity tests

The stationarity of the hydraulic load time series was assessed using the Augmented Dickey-Fuller (ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests, at a significance level of 5%. The ADF test initially indicated non-stationarity of the original series, with a p-value of 0.286026, which is higher than the significance level. After one differencing (d=1), the p-value decreased to 0.129114, still indicating non-stationarity. Only after two differencings (d=2) did the p-value reach 7.300095e-13, confirming the stationarity of the series.

The KPSS test, which assumes stationarity as the null hypothesis, also corroborated the need for differencing. With an initial p-value of 0.100000, the test did not reject the null hypothesis of stationarity, but the combination with the ADF test and the graphical analysis of the twice-differenced series, which showed stable behavior around the mean, confirmed the adequacy of the parameter d=2 for modeling. Stationarity is a prerequisite for applying ARIMA and SARIMA models, ensuring that the mean and variance of the series remain constant over time (Morettin and Toloi, 2018).

Estimation and fitting of the models

The estimation of SARIMA and SARIMAX models was performed using the `auto_arima` function from the `pmdarima` package, seeking to minimize the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), which assess model parsimony and goodness of fit. Initially, SARIMA models, without exogenous variables, were estimated for comparison. The SARIMA(1,2,0)(2,0,1)[12] model showed an AIC of -15,377, while the SARIMA(3,2,0)(2,0,0)[12] had an AIC of -18,023, indicating a better fit.

The introduction of exogenous variables in SARIMAX models demonstrated a significant improvement in performance. The SARIMAX(0,2,3)(2,0,0)[12] model, with all exogenous variables from Table 1, achieved an AIC of -38,882, a substantial drop compared to pure SARIMA models. However, this model presented autocorrelation in the residuals, as indicated by the Ljung-Box test (p<0.05), suggesting that not all information was adequately captured.

The analysis of the coefficients of the SARIMAX models with all exogenous variables revealed that some SPI variables and `pluv_acum2in`, as well as certain autoregressive and moving-average terms, were not statistically significant. This observation can be attributed to multicollinearity among the predictor variables, a common phenomenon in time series (Fávero and Belfiore, 2024). Models with exogenous variables such as `pluv_acum3` (accumulated precipitation in the previous 3 months) and `spi1` (1-month standardized precipitation index) showed the lowest errors, which is consistent with the aquifer’s lagged response to precipitation.

The SARIMAX(2,2,3)(2,0,1)[12] model with the `pluv_acum3` variable stood out with a Mean Absolute Error (MAE) of 0.1041 meters for a 12-month forecast, considered satisfactory for the study’s objective. Another notable model was SARIMAX(2,2,2)(2,0,0)[12] with `spi1`, which achieved an AIC of -52.144 and BIC of -31.066, indicating a robust fit. The final selection of models for fitting and residual diagnostics was based on the AIC, BIC criteria, and the Ljung-Box statistic, prioritizing parsimony and the statistical significance of the coefficients.

Diagnosis of residuals

The diagnosis of the residuals of the adjusted models was performed to verify their representativeness, as recommended by Morettin and Toloi (2018). The evaluation of the scatter plots of the standardized residuals of all adjusted models revealed a random distribution around the zero mean over time, indicating homoscedasticity. This characteristic was formally confirmed by the Heteroskedasticity statistic, which, according to the documentation of the `Statsmodels` library, is analogous to the Goldfeld-Quandt test, whose null hypothesis is the non-heteroscedasticity of the residuals.

However, the Quantile-Quantile plots indicated that the distribution of residuals across all models did not fully adhere to normality, possibly due to the presence of tail values. This non-adherence was corroborated by the Jarque-Bera and Shapiro-Wilk statistical tests, which presented p-values below 0.05 for all adjusted models. Although the non-normality of residuals does not invalidate the models, it may have implications for the determination of confidence intervals (Fávero and Belfiore, 2024). The transformation of the dependent variable, such as the Box-Cox transformation, could mitigate this problem, but it was not explored in the present work.

In summary, the results demonstrated the applicability of statistical machine learning models in predicting the hydraulic load in a contaminated area, with emphasis on SARIMAX models that incorporate exogenous precipitation variables. The ability to predict the water level with a mean absolute error of 0.1041 meters, using accumulated precipitation over 3 months, represents a robust proof of concept. This quantitative evidence is crucial for supporting remediation planning and long-term management of hydrocarbon-contaminated areas, especially in contexts of LNAPL trapping and seasonal water level fluctuations.

4. Conclusion

This work demonstrated the applicability of machine learning models in predicting hydraulic head time series in an area contaminated by aviation kerosene, providing a basis for decision-making regarding the continuation of remediation by hydraulic techniques. It was observed that the aquifer’s hydraulic head exhibits clear seasonal patterns and a lagged response to rainfall. The estimation of SARIMAX models, which incorporated exogenous precipitation variables, significantly improved predictive capacity compared to pure SARIMA models. A SARIMAX model, using three-month accumulated precipitation, achieved a mean absolute error of 0.1041 meters for a twelve-month forecast, which is considered satisfactory. This result established a robust proof of concept, offering a crucial quantitative line of evidence for remediation planning and long-term management of hydrocarbon-contaminated sites, especially in scenarios of LNAPL trapping and seasonal water table fluctuations.

Despite the relevance of the findings, the study presented limitations, such as the linear interpolation of 9.40% of the missing data, whose impact on the models’ sensitivity was not evaluated, and the lack of full adherence of the residuals to normality, which may influence the determination of confidence intervals. Multicollinearity was also observed among some exogenous SPI variables, resulting in non-significant coefficients for certain variables. For future studies, it is suggested to deepen the investigation into multicollinearity and the use of clusters identified by the DBScan algorithm in model estimation. It is also recommended to test approaches such as the Global Forecasting Method and explore more advanced models, such as LSTM, NeuralProphet, and those based on Transformers, to capture non-linearities and improve predictive capacity in specific piezometric wells.

Bibliographic References

Alden, D.F.; García-Rincón, J.; Rivett, M.O.; Wealthall, G.P.; Thomson, N.R. 2024. Complexities of Petroleum Hydrocarbon Contaminated Sites. p. 1-19. In: García-Rincón, J.; Gatsios, E.; Lenhard, R.J.; Atekwana, E.A.; Naidu, R. Advances in the Characterisation and Remediation of Sites Contaminated with Petroleum Hydrocarbons. Environmental Contamination Remediation and Management. Springer, Cham.

Companhia Ambiental do Estado de São Paulo [CET

Companhia Ambiental do Estado de São Paulo [CETESB], 2025. Sistema de Áreas Contaminadas e Reabilitadas. Relatório de Áreas Contaminadas e Reabilitadas no Estado de São Paulo. Disponível em: < https://sigam.ambiente.sp.gov.br/sigam3/Default.aspx?idPagina=17653>. Acesso em: 10 out. 2025.

Organização das Nações Unidas para a Educação, a Ciência e a Cultura [UNESCO]. 2022. UN World Water Development Report 2022: Groundwater: Making the invisible visible. Disponível em: <https://www.unwater.org/publications/un-world-water-development-report-2022>. Acesso em: 12 out. 2025.

Article originating from the Final Course Work of the Specialization in Data Science and Analytics of the MBA USP/Esalq

To learn more about the course, click here and access the MBX Academy platform

You may also like

Economy, Investments And Banking

October 23, 2026

Financial Education, Behaviors and Biases in Financial Decision Making

Personal financial management represents a growing challenge, influenced by cognitive biases that affect decision-making, even in individuals with high knowledge. The present study analyzed the relationship between the level of financial education and the occurrence of behavioral biases in the financial decision-making of adults residing in Brazil. An applied, quantitative, and descriptive-correlational field survey was conducted by administering a structured questionnaire to 163 respondents. The instrument included sociodemographic variables, financial knowledge and habits, and hypothetical situations to identify cognitive biases. The results indicated that higher levels of financial education were associated with better financial control, saving, and investment practices, with significant correlations with income and schooling. However, the presence of behavioral biases (anchoring, loss aversion, scarcity, optimism/negativity, present bias, status quo, overconfidence, and mental accounting) persisted even among individuals with higher schooling and income, with no statistically significant association with the analyzed sociodemographic variables. It was concluded that financial education, although essential for financial organization and control, was not sufficient to neutralize the influence of cognitive biases. The incorporation of behavioral approaches into financial education programs is suggested.

Keywords: Financial behavior; Behavioral economics; Personal finance; Cognitive heuristics; Cognitive biases.

October 22, 2026

Consumer Journey in the Decision to Purchase Infoproducts in the Digital Environment

The expansion of e-commerce and technological transformation have driven the need to understand consumer behavior regarding infoproducts, whose perceived experience has become crucial for professional development. The study aimed to analyze the factors that influenced the purchase decision of infoproducts in the digital environment, considering the stages of the consumer journey. An applied, quantitative, and descriptive research was conducted through a field survey with 100 consumers from Serra Gaúcha. Data collection occurred via a structured questionnaire, using non-probabilistic sampling. The results revealed a qualified audience with strong “mobile-first” behavior. It was found that the availability of detailed information and social approval were the main attraction triggers in the pre-purchase phase. In the transactional and subsequent stages, the minimization of operational friction, the agility of immediate access, and system security proved crucial for conversion and satisfaction, although gaps in human support efficiency were identified. Retention was based on the trust generated by positive past experiences, overcoming the price factor. It was concluded that the infoproduct consumer journey operates as an integrated system, requiring alignment between technological infrastructure, content curation, and support humanization to ensure engagement.

Keywords: E-commerce; Consumer behavior; Consumer experience; Loyalty; Virtual business.

October 22, 2026

Forecast of Aggregate Default in the Brazilian Financial System by Sarimax Model with Macroeconomic Variables

Aggregate non-payment in the Brazilian financial system is a crucial indicator of credit risk, the prediction of which is regulated by CMN Resolutions n.° 4.557/2017 and 4.966/2021 for the prospective measurement of expected losses. Aggregate non-payment was modeled and predicted using the SARIMAX model, projecting its trajectory until December 2026 under optimistic, base, and pessimistic scenarios. 164 monthly observations (July/2012 to February/2026) and macroeconomic variables such as unemployment, Selic, IPCA, IBC-Br, and household debt were used. The specification with contemporaneous variables was selected based on information criteria and coefficient stability, and residual diagnostics confirmed statistical adequacy, with a MAPE of 3.12%. Historical non-payment ranged between 4.33% and 7.64%. The estimation identified unemployment (β = 0.2534; p < 0.001) and Selic (β = 0.2128; p < 0.001) as the main positive determinants, with strong temporal persistence (AR(1) = 0.8965). Selic showed increasing correlation with non-payment (r = 0.62 at lag 12), while unemployment demonstrated a contemporaneous effect (r = 0.33). Projections indicated non-payment between 6.30% and 7.98% in December 2026, with a range of 1.68 percentage points, highlighting direct implications for provisioning in financial institutions.

Keywords: Unemployment; Predictive modeling; Monetary policy; Loss provisioning; Time series.