Mapping and modelling of the spatial distribution of geochemical anomalies is a key task for geochemical exploration. This case study explores mapping and modelling of Cu and Fe geochemical anomalies combining the Sequential Gaussian Co-simulation (SGCS) method and Local Singularity Analysis (LSA) in the Malatya region, SE Turkey. A total of 652 topsoil samples were collected on a regular grid design from the study area. A cell declustering technique was applied to the raw data. Both geometric and zonal anisotropy exist in the directional variogram of Cu and Fe. SGCS was applied to create maps representing an equally probable spatial distribution of Cu and Fe geochemical anomalies. SGCS results have been validated by a variety of tests including the reproduction of the variograms, histograms, descriptive statistics and contour plots. LSA is based on a sliding window estimation approach applied to the grid data created for SGCS realization. LSA was carried out and mapped via the results of the realization on the same grid layout. Critical anomaly threshold values of the variables were identified using the singularity and quantile plot. Integration of SGCS and LSA results provides quantification of the uncertainty of spatial distribution and determination of Cu and Fe critical thresholds.