Joint inversion of flow and seismic data for reservoir parameters is a challenging task in that these disparate datasets are sensitive to different physics and model resolutions for the forward problem. The inverse problem is highly non-linear introducing additional complexity. To overcome some of these challenges we have developed a global optimization method based on very fast simulated annealing (VFSA) and a pilot point based model parameterization scheme. Reservoir simulation is used to create the saturation and pressure distribution with time. The simulation results, are converted to seismic properties using an appropriate rock physics model. Seismic modeling is used to create the seismic response. The objective function is defined as a weighted sum of data misfit and prior model misfit and VFSA is used to derive optimal model parameters. Our results from synthetic examples reveal that the VFSA optimization scheme is robust and pilot point model parameterization is able to obtain reasonable descriptions of the reservoir. We further propose a probability based pilot point parameterization, where prior knowledge is used to compute the probability to draw the pilot points. In this way, the model parameters can be reduced further. To incorporate the small scale heterogeneity, we combine the pilot point based inversion method with sequential Gaussian simulation to create stochastic models.