We present a model of magma bodies formation by injection of hot rhyolitic magmatic dikes leading to their incremental accumulation into the plutons and magma chambers in the upper and middle crust. Our 2D model simulates random or organized dike injection into a selected rock volume, calculation of magma and rock displacement based on analytical solution of an elastic problem of elliptical cavity expansion, a realistic melt phase diagrams for country rock and magma. Lagrangian particle transport is calculated in order to reduce numerical dissipation and avoid unphysical mixing. Thermal histories in individual batches of magma and country rocks are recorded. We further combine this model with Bindeman and Melnik (2016) zircon crystallization/dissolution software and compute zircon survival histories in each Lagrangian particle. The model predicts shapes of realistic T-t histories, zircon age distribution in different portions within a progressively growing reservoir and generates output to estimate crustal vs mantle contributions (e.g. Hf and O isotopes in zircons).Simulations reveal that the rate of melt production is highly variable in space and time, eruptible magma batches form in clusters, period of initial magmatic incubation is followed by crustal rock melting and formation of a large volume of eruptible magma with high melt fraction. Zircon survival and host rock eruptibility depends on magma injection duration. After 700 y only ~2 vol % of molten rock can be erupted, but due to slow dissolution most of zircon crystal are preserved. After 1500 years eruptible rock amount reaches 8% but significant number of zircons looses their age information.