Mathematical programming has been widely applied for the planning of natural gas production infrastructure development. As the production infrastructure involves large investments and is expected to remain in operation over several decades, the factors that will impact the gas production but cannot be foreseen before the development of the infrastructure need to be taken into account in the planning. Therefore, two scenario-based two-stage stochastic programming models are developed to facilitate natural gas production infrastructure development under uncertainty. One is called the stochastic pooling model, which tracks the qualities of gas streams throughout the production network via a generalized pooling model. The other is an enhancement of the stochastic pooling model with the consideration of pressure. Either model results in a large-scale nonconvex mixed-integer nonlinear programming (MINLP) problem, for which a global optimal solution, although very important for a problem that involves large investments, is very difficult to obtain. A novel optimization method, called nonconvex generalized Benders decomposition (NGBD), is developed for efficient global optimization of the large-scale nonconvex MINLP. Case studies of a real industrial natural gas production system show the advantages of the proposed stochastic programming models over deterministic optimization models, as well as the dramatic computational advantages of NGBD over a state-of-the-art global optimization solver.