Processed-based dynamic vegetation models (DVMs) are important component of earth system modelling and important for assessing future climate impacts. The DVMs are often highly complex and non-linear with many tuning parameter. By matching model output to observed CO2 fluxes, we aim to estimate the parameters. Stochastic optimization methods have proved useful for similar optimization problems. These methods use Monte Carlo integration to propagate a probabilistic model that localise regions of high quality solutions. The integration requires multiple evaluation of the DVM, which are very expensive to compute. To reduce the number of cost evaluation we propose an improved optimization method that uses previous samples to improve convergence and reduce the number of function evaluations.