Approaches to lake water quality modeling are critically examined with particular attention to the formulation of water quality transport as the link between hydrodynamics and biogeochemical reaction. A linked water quality model for shallow lakes includes three major components: a biogeochemical reaction component, a lake hydrodynamics component and a water quality transport component. State-of-the-art modeling approaches for each component are reviewed, and a synopsis of phosphorus dynamics in shallow lakes is given. For the water quality transport component, review of the literature shows two significantly different approaches to water transport: a lumped component approach based upon multiple fully-mixed boxes, and a continuum approach employing the finite difference method to approximate the continuous governing equations. The multiple-box model is shown in an analysis of the kindred fully-mixed tanks-in-series conceptual reactor model to create an excessive implicit dispersion due to its formulation. This leads to a model in which the model mass transport is not closely related to the properties of the physical system. being modeled. Rather, dispersive transport in the model is shown to depend heavily upon the model formulation -- the model transport parameters thus cannot be specified from hydrodynamic data but must be found by calibration. In direct contrast, the finite difference model maintains a far closer approximation to the physical system and permits direct specification of water quality transport from the actual lake hydrodynamics. To support these arguments, a computer program incorporating both a multiple-box model and an alternative one-dimensional finite difference model is developed and applied to Lake Balaton in Hungary. The biogeochemical component of both models is a four component phosphorus-phytoplankton interaction model originally proposed by van Straten (1980). Hydrodynamic information is supplied to the one-dimensional finite difference model by linkage to a transient two-dimensional single-layer model of wind-driven circulation. A one-dimensional dispersion coefficient is computed from the two-dimensional velocity field using a method based upon that of Fischer (1966, 1969) and Holley, Harleman and Fischer (1970), but proceeding from the assumption that advection rather than turbulent diffusion dominates lateral mixing. The finite difference model developed for Lake Balaton consists of forty grids and is used to simulate a representative period from early spring to late summer. The results are contrasted with those produced by a four-box model of the lake using long-term average advection and calibrated dispersive exchange flows. The finite difference model is found to lead to a predicted behavior more similar to that observed in field data collected from Lake Balaton. Experimentation with the models is conducted to examine the behavior of the lake and the dominant factors leading to that behavior.