This paper presents a novel concept in the thermodynamic derivation of phase diagrams for clay minerals that incorporates fuzzy transition zones. This new technique yields phase diagrams that have graded (fuzzy) zones of mineral occurrences and includes compositional variability within mineral groups. For the construction of these diagrams, 170 minerals belonging to nine different subgroups were used, based on a fuzzy mathematical description of their ‘grades’ or ‘belonging-ness’. Standard free energies of formation of all the minerals were derived and all possible pairs of mineral equilibria were evaluated. Relative intensities of mineral occurrences were determined and membership values of each type of mineral at various zones in a 2D or 3D space were graphically represented. Computations and graphical representations were carried out using programs developed in Mathematica. Diagrams were derived for 25°C, 1 bar with a solution phase containing Si(OH)4, K+, Na+, H+, Ca2+ and Mg2+ under conditions of gibbsite, goethite and ferrous oxide saturation. The resulting diagrams, unlike conventional phase diagrams, show multimineral assemblages, with varying occurrences of different minerals and provide a realistic representation of clay mineral occurrences formed by surface geochemical processes. They show that on the Earth’s surface, only montmorillonite can almost completely predominate the inorganic phase followed by kaolinite, illite and beidellite. Nontronite, glauconite, celadonite and vermiculite would not be neoformed in substantial amounts. A general conformity of derived phase equilibria with experimentally observed equilibria is also observed.