This work investigates upscaling of capillary pressure curves for modeling gravity segregation under the influence of capillary heterogeneity. We consider two flowing phases driven by gravity and capillary forces and seek the saturation spatial and temporal variation until equilibrium is reached. Existing upscaling methods, found in the literature, are applied to a number of cases. Resulting saturation solutions are compared to fine-scale simulations and the different methods are evaluated, showing in general that the popular capillary limit method produces the best results. However, a large number of cases are found to have significant errors. This leads to the conclusion that existing methods are often inadequate. We therefore propose a new optimization-based upscaling method. Using this approach it is shown that capillary pressure can be upscaled to Brooks-Corey type functions and produce accurate upscaled simulations matching the fine-scale solutions. Optimization upscaling is computationally demanding and requires the fine-scale simulations for minimizing an objective function. However, it is shown that the upscaled curves are robust, i.e., they can be applied to different permeability realizations to calculate ensemble average saturation solutions.