One of the major problems in the definition of early warning systems at active volcanoes is represented by the need of understanding the relationships between changes in quantities measured at the Earth surface, and movement of magma at depth. Among the measured quantities, surface deformation and gravity are being measured with increasingly frequency and through innovative methods which include ground-based and remote sensing techniques. In order to contribute to establishing a link between measurements and on-going processes at active volcanoes, we have developed a forward approach which consists in modelling the time-space-dependent dynamics of magma flow in the sub-surface volcanic environment (magma reservoirs and volcanic conduits or dykes) and associated ground deformation, seismic signals, and gravity changes. Such an approach provides a consistent view of magma and rock dynamics associated with a variety of processes occurring below the volcano and potentially capable to generate a volcanic eruption, such as magma chamber replenishment, mixing of compositionally different magmas, volatile exsolution, magma ascent in a volcanic conduit or dyke, etc. The numerical code for magma flow (Longo et al., Geophys. Res. Lett. 2006) describes the time-dependent 2D dynamics of a compressible-to-incompressible homogeneous multicomponent mixture made of liquid in equilibrium with an H2O+CO2 gas phase at local P-T-X conditions. The finite element numerical algorithm consists in a time-space discretization with Galerkin least-squares and discontinuity capturing terms, which allow high numerical stability. The code is implemented with the most updated models for non-ideal multicomponent gas-liquid thermodynamics and for the relevant properties density and viscosity, under the assumption of Newtonian magma rheology. The numerical code for rock dynamics (O’Brien and Bean, Geophys. Res. Lett. 2004) is based on an elastic lattice method which consists of particles arranged on a cubic lattice which interact through a central force term and a bond-bending force, and which allows taking into account the topographic surface and rock heterogeneities. One-way coupling between the two models occurs through the definition of the time-space force source for rock movements defined by the time-space stress distribution at the magma-rock boundary calculated by the fluid dynamic model. We present here the results of a study on magma chamber replenishment in a possible shallow chamber at Campi Flegrei. The system initial and boundary conditions are defined on the basis of the whole investigation carried out in the frame of a 2-year INGV-Department of Civil Protection project which included several tens of geophysicists, petrologists, geochemists, field volcanologists, and numerical modellers. The system is schematized as a 2D long dyke feeding a shallow and small magma chamber. While the chamber initially hosts a magma with phonolitic composition, the dyke provides CO2-rich magma with trachytic composition. The simulations show the dynamics of light magma plume rise, convection, and mixing with the resident magma, and other dynamic aspects of such processes which were not known before. Parametric studies allowed us to vary the geometry of the chamber, the volatile contents of the two magmas, and the initial overpressure at dyke base. By integrating the contribution from any point in space in the simulated domain, we model the expected free-air corrected gravity change associated with the simulated dynamics. The computed gravity changes are negative, due to the replacement at shallow depth of the denser phonolitic magma by the CO2-rich, lighter trachytic magma, with maximum values of some tens of mugal above the chamber. The numerical simulations of rock dynamics employ either homogeneous or heterogeneous rock properties, the latter defined on the basis of the results of seismic tomography experiments at Campi Flegrei (Zollo et al., Geophys. Res. Lett. 2003). A subset of simulations also employ a schematized 2D topography of the caldera region. The results show the time-space distribution of ground deformation and the seismic signals obtained by high-pass filtering at 100 s according to the frequency band of typical broadband seismometers. Ground movement shows oscillations related to the dynamics of magma convection inside the chamber, with maximum vertical displacements of the order of cm and typical periods concentrating in two frequency bands corresponding to several tens and a few to a few tens of s. A number of characteristics of the quasi-static and seismic-band ground movements associated with magma chamber replenishment at Campi Flegrei are described. In order to be fully used for the development of early warning systems at active volcanoes, our numerical codes and simulations require the development of more accurate techniques to detect and describe geodetic signal variations over the time-scale of tens of seconds to minutes.