Water quality management under the watershed approach of Total Maximum Daily Load (TMDL) programs requires that water quality standards be maintained throughout the year. The main purpose of this research was to develop a methodology that incorporates inter-temporal variations in stream conditions through statistical distributions of pollution loading variables. This was demonstrated through a cost minimization mixed-integer linear programming (MIP) model that maintains the spatial integrity of the watershed problem. Traditional approaches for addressing variability in stream conditions are unlikely to satisfy the assumptions on which these methodologies are founded or are inadequate in addressing the problem correctly when distributions are not normal. The MIP model solves for the location and the maximum capacity of treatment plants to be built throughout the watershed which will provide the optimal level of treatment throughout the year. The proposed methodology involves estimation of parameters of the distribution of pollution loading variables from simulated data and use of those parameters to re-generate a suitable number of random observations in the optimization process such that the new data preserve the same distribution parameters. The objective of the empirical model was to minimize costs for implementing pH TMDLs for a watershed by determining the level of treatment required to attain water quality standards under stochastic stream conditions. The output of the model was total minimum costs for treatment and selection of the spatial pattern of the least-cost technologies for treatment. To minimize costs, the model utilized a spatial network of streams in the watershed, which provides opportunities for cost-reduction through trading of pollution among sources and/or least-cost treatment. The results were used to estimate the costs attributable to inter-temporal variations and the costs of different settings for the 'margin of safety'. The methodology was tested with water quality data for the Paint Creek watershed in West Virginia. The stochastic model included nine streams in the optimal solution. An estimate of inter-temporal variations in stream conditions was calculated by comparing total costs under the stochastic model and a deterministic version of the stochastic model estimated with mean values of the loading variables. It was observed that the deterministic model underestimates total treatment cost by about 45 percent relative to the 97th percentile stochastic model. Estimates of different margin of safety were calculated by comparing total costs for the 99.9th percentile treatment (instead of an idealistic absolute treatment) with that of the 95th to 99th percentile treatment. The differential costs represent the savings due to the knowledge of the statistical distribution of pollution and an explicit margin of safety. Results indicate that treatment costs are about 7 percent lower when the level of assurance is reduced from 99.9 to 99 percent and 21 percent lower when 95 percent assurance is selected. The application of the methodology, however, is not limited to the estimation of TMDL implementation costs. For example, it could be utilized to estimate costs of anti-degradation policies for water quality management and other watershed management issues.