A multispecies numerical code was developed to simulate flow and mass transport with kinetic adsorption in variable-density flow systems. The two-dimensional code simulated the transport of bromide (Br⁻), a nonreactive tracer, and lithium (Li⁺), a reactive tracer, in a large-scale tracer test performed in a sand-and-gravel aquifer at Cape Cod, Massachusetts. A two-fraction kinetic adsorption model was implemented to simulate the interaction of Li⁺ with the aquifer solids. Initial estimates for some of the transport parameters were obtained from a nonlinear least squares curve-fitting procedure, where the breakthrough curves from column experiments were matched with one-dimensional theoretical models. The numerical code successfully simulated the basic characteristics of the two plumes in the tracer test. At early times the centers of mass of Br⁻ and Li⁺ sank because the two plumes were closely coupled to the density-driven velocity field. At later times the rate of downward movement in the Br⁻ plume due to gravity slowed significantly because of dilution by dispersion. The downward movement of the Li⁺ plume was negligible because the two plumes moved in locally different velocity regimes, where Li⁺ transport was retarded relative to Br⁻. The maximum extent of downward transport of the Li⁺ plume was less than that of the Br⁻ plume. This study also found that at early times the downward movement of a plume created by a three-dimensional source could be much more extensive than the case with a two-dimensional source having the same cross-sectional area. The observed shape of the Br⁻ plume at Cape Cod was simulated by adding two layers with different hydraulic conductivities at shallow depth across the region. The large dispersion and asymmetrical shape of the Li⁺ plume were simulated by including kinetic adsorption-desorption reactions.