A geochemical model was incorporated into a stratification model for lakes to create the model package: DYCD-CORE, a numerical code that couples the thermal and hydrodynamic capabilities of DYRESM and the geochemical capabilities of the reactive transport model CORE2D V4. Based on the chemical composition of solutes calculated in each node for each time step, density was computed using specific partial molal volumes of all considered solutes and fed back into the stratification module of the program package. The density calculated each time step leads to a strong coupling of hydrodynamics and hydrogeochemistry and reflects the complex interaction as it is present in all lakes. To demonstrate the functionality of the numerical approach, an example of an iron-meromictic lake was chosen, where the reactivity of the dissolved iron kept the water body perennially stratified. To critically validate the model results, temperatures were continously measured at high vertical and temporal resolution in a field investigation of Waldsee (near Dobern, Germany). Multiparameterprobe profiles and water samples confirmed the continous chemical stratification and served as initial and boundary conditions for the simulation period. The model package DYCD-CORE could reproduce the permanent stratification as it were in the lake. A demonstration run using the standard UNESCO equation for density, and hence assuming non-reactive solutes, failed entirely. Hence, stratification models using salinity for density are not suited for simulating density created by lake-internal geochemical transformation of solutes. However, density can be based directly on the simultaneous numerical simulation of lake geochemistry. Predictive modeling of changing lake circulation in a variable climate or considering change of use will require a proper inclusion of the geochemistry as demonstrated in this paper.