This study presents a systematic geochemical modeling methodology using a combination of data analysis techniques suitable for the predictive modeling of metal contamination processes in altered, heterogeneous fluvial sediment system as contamination receptor. Unlike conventional studies, this approach first characterizes the receptor fluvial sediments using (1) major ions (Fe, Al, K, Na, Ca, Mg, P), indicative of sediment mineralogy; and (2) standard soil parameters (pH, electrical conductivity, loss-on-ignition, CaCO3 content, clay content), indicative of sediment chemistry, controlling metal behavior (e.g. adsorption, precipitation). In addition to unsupervised hierarchical cluster analysis and Q-mode (between-sample) principal component analysis, multivariate homogeneity tests such as the Mann-Whitney and Kruskal-Wallis median tests identify and verify distinct, geochemically homogeneous areas which, in this study, correspond to the contemporary riverbed sediments, the post-river regulation active alluvial plain, and to the old terrace developed under braided-meandering pre-river regulation conditions of the Drava River. The main geochemical processes in the receptor floodplain soils and sediments, characterized by associations of major ion concentrations and soil parameters, modeled with regression, partial correlation analysis and correlograms are: (1) calcareous geochemistry: Ca-Mg-pH-(CaCO3 content), in the stream and alluvial plain sediments, defined by the regional geochemical background, (2) clay minerals (sericite-illite): K-Al-Fe-(Mg, clay content), in the river terrace soils, while (3) organic matter: loss-on-ignition-electrical conductivity-P, dominates the topsoil geochemistry in the floodplain. Second, the contamination source is modeled by geochemical fingerprinting using (1) metal associations, characteristic to the upstream mines of the Mississippi Valley-type Pb-Zn-Cd deposits such as the minerals wulfenite Pb(MoO4), descloizite PbZn(VO4)(OH), and pyromorphite Pb5(PO4)3Cl, and (2) metal associations, characteristic to mafic-ultramafic rocks in the catchment area and to the regional technogenic (smelting and metallurgy) geochemical signals (Ni-Cr-Co-V-(Fe)). This step uses Enrichment Factor calculation, regression and partial correlation analysis. Third, the metal contamination-receptor sediment interaction is modeled with regression analysis, applied to distinct sediment geochemistry areas, complemented by thermodynamic hydro-geochemical modeling. Finally, a simple but efficient predictive modeling procedure is developed to assess the fate of metal contamination introduced into the heterogeneous fluvial system: stepwise supervised R-mode (between-parameter) cluster analysis is used where metals are added one-by-one to the previously defined sediment-soil parameter clusters (calcareous geochemistry, clay minerals, organic matter). Accordingly, the Pb-Zn-Cd metals associate with carbonates in the stream and alluvial plain sediments, while Ni-Cr-Co-V-(Fe) metals associate with clays, traceable predominantly in the river terrace, just like heavy minerals. Metal(loid)s also tend to be adsorbed to the reactive surfaces of clay minerals and organic matter, especially in river terrace soils. It is concluded that the applied systematic geochemical data analysis and modeling methodology provides an efficient tool for the characterization and modeling of the complex processes of regional historical metal contamination impacting a heavily disturbed soil-sediment system of a large-river fluvial terrain.