In siliciclastic and carbonate reservoirs, depositional facies are commonly described as being organized in cyclic successions overprinted by diagenesis. Most reservoir modeling workflows are not able to stochastically reproduce such patterns. Herein, a novel geostatistical method is developed to model depositional facies architectures that are rhythmic and cyclic, together with superimposed diagenetic facies. The method uses truncated Pluri-Gaussian random functions constrained by transiograms. Cyclicity is defined as an asymmetric ordering between facies, and its direction is given by a three-dimensional vector, called shift. This method is illustrated in two case studies. Outcrop data of the Triassic Latemar carbonate platform, northern Italy, are used to model shallowing-upward facies cycles in the vertical direction. A satellite image of the modern Bermuda platform interior is used to model facies cycles in the windward-to-leeward lateral direction. Because depositional facies architectures are modeled using two Gaussian random functions, a third Gaussian random function is added to model diagenesis. Thereby, depositional and diagenetic facies can exhibit spatial asymmetric relationships. The method is applied in the Latemar carbonate platform that experiences syndepositional dolomite formation. The method can also incorporate proportion curves to model nonstationary facies proportions. This is illustrated in Cretaceous shallow-marine sandstones and mudstones of the Book Cliffs (Utah), for which cyclic fades and diagenetic patterns are constrained by embedded transition probabilities.