In this work, a crystallographic thermal creep model is proposed for Zr alloys that accounts for the hardening contribution of solutes via their time-dependent pinning effect on dislocations. The core-diffusion model proposed by Soare and Curtin (2008a) is coupled with a recently proposed constitutive modeling framework (Wang et al., 2017, 2016) accounting for the heterogeneous distribution of internal stresses within grains. The Coble creep mechanism is also included. This model is, in turn, embedded in the effective medium crystallographic VPSC framework and used to predict creep strain evolution of polycrystals under different temperature and stress conditions. The simulation results reproduce the experimental creep data for Zircaloy-4 and the transition between the low (n~1), intermediate (n~4) and high (n~9) power law creep regimes. This is achieved through the dependence on local aging time of the solute-dislocation binding energy. The anomalies in strain rate sensitivity (SRS) are discussed in terms of core-diffusion effects on dislocation junction strength. The mechanism-based model captures the primary and secondary creep regimes results reported by Kombaiah and Murty (2015a, 2015b) for a comprehensive set of testing conditions covering the 500 to 600 o C interval, stresses spanning 14 to 156 MPa, and steady state creep rates varying between 1.5• 10-9 s-1 to 2• 10-3 s-1. There are two major advantages to this model with respect to more empirical ones used as constitutive laws for describing thermal creep of cladding: 1) specific dependences on the nature of solutes and their concentrations are explicitly accounted for; 2) accident conditions in reactors, such as RIA and LOCA, usually take place in short times, and deformation takes place in the primary, not the steady-state creep stage. As a consequence, a model that accounts for the evolution with time of microstructure is more reliable for this kind of simulation.