Option pricing based on stochastic volatility models can reproduce volatility skew/smile observed on the market. The use of stochastic volatility models typically relies on either closed-form solutions or Monte Carlo simulations, each with obvious limitations. Finite difference methods to solve the resulting time-dependent PDE in two space dimensions provide a powerful alternative, with advantages such as natural accommodation of the correlation and variable coefficients. One less obvious but potentially critical issue is the implementation of boundary condition imposed at the zero volatility boundary. In this work, we consider the full PDE problem and associate different boundary conditions with the corresponding stochastic processes in various regimes, and use high-order finite difference methods to analyze the boundary impact. Using the settings of Heston and SABR models, we present numerical results, in terms of the market observable "volatility smile" curve, to demonstrate the ramifications of the boundary treatments. The importance of this issue is made obvious in exotic options such as the one-touch options in currency trading. Comparisons with other approaches such as Monte Carlo simulations are also performed to demonstrate the particular features of the finite difference methods. (Joint with Peter Laurence)