The subsurface transport of phosphorus introduced by the disposal of treated sewage effluent to ground-infiltration disposal beds at the Massachusetts Military Reservation on western Cape Cod was simulated with a three-dimensional reactive-transport model. The simulations were used to estimate the load of phosphorus transported to Ashumet Pond during operation of the sewage-treatment plant from 1936 to 1995 and for 60 years following cessation of sewage disposal. The model accounted for spatial and temporal changes in water discharge from the sewage-treatment plant, ground-water flow, transport of associated chemical constituents, and a set of chemical reactions, including phosphorus sorption on aquifer materials, dissolution and precipitation of iron- and manganese-oxyhydroxide and iron phosphate minerals, organic carbon sorption and decomposition, cation sorption, and irreversible denitrification. The flow and transport in the aquifer were simulated by using parameters consistent with those used in previous flow models of this area of Cape Cod, except that numerical dispersion was much larger than the physical dispersion estimated in previous studies. Sorption parameters were fit to data derived from phosphorus sorption and desorption laboratory column experiments. Rates of organic carbon decomposition were adjusted to match the location of iron concentrations in an anoxic iron zone within the sewage plume. The sensitivity of the simulated load of phosphorus transported to Ashumet Pond was calculated for a variety of processes and input parameters. Model limitations included large uncertainties associated with the loading of the sewage beds, the flow system, and the chemistry and sorption characteristics in the aquifer. The results of current model simulations indicate a small load of phosphorus transported to Ashumet Pond during 1965-85, but this small load was particularly sensitive to model parameters that specify flow conditions and the chemical process by which non-desorbable phosphorus is incorporated in the sediments. The uncertainties were large enough to make it difficult to determine whether loads of phosphorus transported to Ashumet Pond in the 1990s were greater or less than loads during the previous two decades. The model simulations indicate substantial discharge of phosphorus to Ashumet Pond after about 1965. After the period 2000-10 the simulations indicate that the load of phosphorus transported to Ashumet Pond decreases continuously, but the load of phosphorus remains substantial for many decades. The current simulations indicate a peak in phosphorus discharge to Ashumet Pond of about 1,000 kilograms per year during the 1990s; however, comparisons of simulated phosphorus concentrations with measured concentrations in 1993 indicate that the peak in phosphorus load transported to Ashumet Pond may be larger and moving more quickly in the model simulations than in the aquifer. The results of the three-dimensional reactive-transport simulations are consistent with the loading history, experimental laboratory data, and field measurements. The results of the simulations adequately reproduce the spatial distribution of phosphorus concentrations measured in 1993, the magnitude of changes in phosphorus concentration with time in a profile near the disposal beds following cessation of sewage disposal, the observed iron zone in the sewage plume, the approximate flow of treated sewage effluent into Ashumet Valley, and laboratory-column data for phosphorus sorption and desorption.