Materials and methods

Multi-source data collection and feature selection

To develop a robust predictive model for the log10-transformed sediment–water pseudo-partitioning coefficient (log_K_d) of ECs, we compiled a multi-source database comprising peer-reviewed publications retrieved via Google Scholar and China National Knowledge Infrastructure, national and regional monitoring networks, and publicly accessible environmental and sediment datasets (Supplementary Table S1). The dataset covered three representative EC classes (PFASs, ABs, and EDCs), comprising 5085 paired sediment–water concentration records from 1093 sampling sites across China's seven major river basins. For each compound at each site, _K_d (L kg−1) was calculated as [29]:where _C_s and _C_w represent the measured EC concentrations in sediment (ng kg−1) and water phase (ng L−1), respectively [29].

Kd=CsCw

To characterize the spatiotemporal variability of _K_d, we incorporated 22 explanatory variables into the ML models and categorized them into three scales. Microscale variables included molecular descriptors from ChemSpider and Chemicalize [30]: molecular weight, number of rotatable bonds, number of hydrogen bond acceptors, maximum molecular projection diameter (_D_max), and log_K_ow. Mesoscale variables included water-quality parameters (i.e., temperature, dissolved oxygen (DO), pH, chemical oxygen demand (COD), NH3-N) and sediment properties (i.e., total organic carbon content and the relative proportions of clay, silt, and sand). These variables were included because organic-matter characteristics, redox conditions, and nutrient-related biogeochemical processes may influence mobility, sorption behavior, and sediment–water partitioning of contaminants [31,32]. Water-quality data were obtained preferentially from in situ measurements at sampling sites or, when unavailable, from the nearest national monitoring stations in the same river basin as proxy environmental background conditions [6]. Although this approach may not fully capture fine-scale spatial heterogeneity within river reaches, it provides a feasible approximation for basin-scale spatiotemporal analysis. Macroscale variables included runoff, precipitation from the National Tibetan Plateau Scientific Data Center, urban wastewater discharge from the China Urban Construction Statistical Yearbook, and land-use patterns from satellite remote sensing imagery [33].

We used a standardized preprocessing workflow to harmonize heterogeneous datasets [34]. EC concentrations were spatiotemporally aligned with environmental and anthropogenic drivers using nearest-neighbor matching based on latitude–longitude coordinates and sampling dates [35]. Numeric variables were normalized using z-score transformation, while categorical land-use variables were quantified as proportional coverage within multi-scale buffers derived from time-series satellite imagery [36,37]. Buffer sensitivity analysis identified a 500 m radius as yielding optimal model performance (Supplementary Fig. S1), which was adopted for subsequent modeling. For sites lacking direct sediment property measurements, values from spatially proximal sites within local basins were applied and validated using clustering consistency analysis (Supplementary Fig. S2) [38]. Before model development, samples with irretrievable missing values in key explanatory variables or with EC concentrations below the analytical detection limits were excluded, using consistent screening criteria to ensure data completeness and model reliability. This screening reduced the sample size but improved consistency and completeness of multi-source environmental variables collected from heterogeneous literature and monitoring datasets. The final dataset comprised 2377 complete paired sediment–water records for model development.

MB-MHA learning framework versus traditional ANN for predicting logKd

We developed a cross-scale MB architecture to capture the hierarchical drivers of EC partitioning in sediment–water systems by explicitly assigning micro-, meso-, and macroscale predictors to independent branches. Each branch comprised two to three dense layers with nonlinear activation functions and dropout regularization, ensuring both representational capacity and prevention of overfitting [26]. Branch-specific outputs were subsequently integrated via a MHA mechanism, which adaptively weighted their relative contributions and facilitated the learning of cross-scale interactions among predictors through a learnable matrix [39]. This design preserved scale-specific information and enabled synergistic integration across different environmental drivers. For benchmarking, baseline ANN models were constructed by concatenating all predictors into a single input stream.

Both the baseline ANN and MB-MHA models were trained using the Adam optimizer, with mean squared error as the loss function. To ensure spatial independence and eliminate site-level data leakage, all observations from a given monitoring location were assigned exclusively to either the training or the test subset [24]. Under this constraint, the entire dataset was randomly split into training and test subsets at an 80:20 ratio. Fivefold cross-validation and model hyperparameter tuning, including learning rate, batch size, and dropout rate, were performed only on the training set using Bayesian optimization implemented with Hyperopt [40]. Predictive performance was evaluated on the independent test set using the coefficient of determination (_R_2) and root-mean-square error [30].

We conducted ablation experiments to quantify the contribution of different feature groups to the model's performance [41]. Feature-level interpretability was further examined using permutation importance analysis, and partial dependence plots (PDPs) were used to illustrate the nonlinear responses and threshold behaviors of key variables on predicted log_K_d values [30]. Branch-wise attention matrices were visualized to examine cross-scale interactions within the MB-MHA architecture [39]. In addition, t-distributed stochastic neighbor embedding (t-SNE) was employed to compare feature embeddings before and after the introduction of the MHA module, revealing enhanced clustering and discriminative structures that support the effectiveness of the MB-MHA framework [42].

Molecular dynamics simulation for mechanistic insights

Perfluorooctanoic acid (PFOA) and perfluorobutanoic acid (PFBA), tetracycline (TC) and ofloxacin (OFX), bisphenol A (BPA) and 4-nonylphenol (4-NP) were selected as representative compounds of PFASs, ABs, and EDCs, respectively, based on their widespread environmental occurrence, high detection frequencies, and potential human-health risks [43]. These selected compounds exhibited distinct physicochemical properties, enabling comparison of sediment–water partitioning behaviors across multiple EC subclasses under different environmental conditions. Their three-dimensional structures (Supplementary Table S2) were converted into GROMACS-compatible topology files using PRODRG, with atom types, charges, and bond parameters refined according to the CHARMM36 force field [27].

The montmorillonite (MMT) unit cell was obtained from the American Mineralogist Crystal Structure database and neutralized with Na+ ions [44]. An Na-MMT interlayer system was constructed by splitting the clay platelet and inserting an interlayer gallery, in which ten EC molecules were randomly positioned, and the remaining pore volume was filled with water. An organo-MMT system was further generated by incorporating humic substances into the interlayer region.

The ClayFF force field was applied to parameterize the MMT, while CHARMM36, which is compatible with ClayFF, was used for humic substances and ECs [45]. Water molecules were modeled using the simple point charge scheme. All simulations were conducted in GROMACS2023. After steepest-descent energy minimization, the systems underwent NVT equilibration followed by a 50 ns production run at 298 K. Temperature was controlled using the V-rescale thermostat, and short-range van der Waals and electrostatic interactions were truncated at 1.5 nm [46]. Long-range electrostatics were treated using the Particle Mesh Ewald method. All simulations employed three-dimensional periodic boundaries with a 1.0 fs time step.

Spatiotemporal forecasting and mapping of logKd for ECs in the Greater Bay Area, China

The Greater Bay Area in southern China (21°30′–24°40′ N, 111°21′–114°53′ E) is one of the world's most densely urbanized and industrialized regions. It contains highly interconnected river networks influenced by intensive anthropogenic emissions and complex hydrological regulation [47]. The Greater Bay Area therefore provides a suitable demonstration region for evaluating spatiotemporal variation in the sediment–water partitioning of ECs.

We used the optimized MB-MHA model to generate monthly log_K_d projections for representative ECs, including PFOA, PFBA, TC, OFX, BPA, and 4-NP, at representative months (January, April, July, and October) in 2026. Time-varying input variables, including conventional water quality parameters (DO, pH, COD, NH3-N, temperature) and meteorological features (runoff and precipitation), were forecasted using long short-term memory (LSTM) networks [24]. Urban wastewater discharge and land-use patterns were predicted using linear regression and least-squares trend fitting, respectively. Sediment properties were treated as temporally invariant due to their limited short-term variability and slow geochemical evolution.

Forecasted drivers were spatially aligned with hydrological units and fed into the MB-MHA model to generate basin-scale log_K_d surfaces at approximately 50 m spatial resolution. The resulting monthly maps capture predicted transitions between high-_K_d zones, where ECs are more likely to accumulate in sediments, and low-_K_d zones, where ECs are more likely to remain mobile in the aqueous phase. These maps provide a quantitative basis for identifying potential contamination hotspots and mobility-driven risk areas in the Greater Bay Area. The regional-scale distribution patterns were further interpreted alongside the micro- and mesoscale controlling factors identified by the ML and MD analyses.