A chance‐constrained stochastic programming model is developed for water quality optimization. It determines the least cost allocation of waste treatment plant biochemical oxygen demand (BOD) removal efficiencies, subject to probabilistic restrictions on maximum allowable instream dissolved oxygen deficit. The new model extends well beyond traditional approaches that assume streamflow is the sole random variable. In addition to streamflow, other random variables in the model are initial in‐stream BOD level and dissolved oxygen (DO) deficit; waste outfall flow rates, BOD levels and DO deficits; deoxygenation k1, reaeration k2, and sedimentation‐scour rate k3 coefficients of the Streeter‐Phelps DO sag model; photosynthetic input‐benthic depletion rates Ai, and nonpoint source BOD input rate Pi for the Camp‐Dobbins extensions to the Streeter‐Phelps model. These random variables appear in more highly aggregated terms which in turn form part of the probabilistic constraints of the water quality optimization model. Stochastic simulation procedures for estimating the probability density functions and covariances of these aggregated terms are discussed. A new chance‐constrained programming variant, imbedded chance constraints, is presented along with an example application. In effect, this method imbeds a chance constraint within a chance constraint in a manner which is loosely associated with the distribution‐free method of chance‐constrained programming. It permits the selection of nonexpected value realizations of the mean and variance estimates employed in the deterministic equivalents of traditional chance‐constrained models. As well, it provides a convenient mechanism for generating constraint probability response surfaces. A joint chance‐constrained formulation is also presented which illustrates the possibility for prescription of an overall system reliability level, rather than reach‐by‐reach reliability assignment.