Assimilation of in situ and satellite data in mechanistic terrestrial ecosystem models helps to constrain critical model parameters and reduce uncertainties in the simulated energy, water and carbon fluxes. So far the assimilation of eddy covariance measurements from flux-tower sites has been conducted mostly for individual sites ("single-site" optimization). Here we develop a variational data assimilation system to optimize 21 parameters of the ORCHIDEE biogeochemical model, using net CO2 flux (NEE) and latent heat flux (LE) measurements from 12 temperate deciduous broadleaf forest sites. We assess the potential of the model to simulate, with a single set of inverted parameters, the carbon and water fluxes at these 12 sites. We compare the fluxes obtained from this "multi-site" (MS) optimization to those of the prior model, and of the "single-site" (SS) optimizations. The model-data fit analysis shows that the MS approach decreases the daily root-mean-square difference (RMS) to observed data by 22%, which is close to the SS optimizations (25% on average). We also show that the MS approach distinctively improves the simulation of the ecosystem respiration (Reco), and to a lesser extent the gross primary productivity (GPP), although we only assimilated net CO2 flux. A process-oriented parameter analysis indicates that the MS inversion system finds a unique combination of parameters which is not the simple average of the different SS sets of parameters. Finally, in an attempt to validate the optimized model against independent data, we observe that global-scale simulations with MS optimized parameters show an enhanced phase agreement between modeled leaf area index (LAI) and satellite-based observations of normalized difference vegetation index (NDVI).