As noted by Yeh and Tripathi [1989], Rathfelder et al. [1991], Reeves and Abriola [1994], and others, nonequilibrium significantly alters model formulation. With LEA, the phase concentrations that govern component fluxes are determined algebraically, locally in space, by component total concentrations. With kinetics, phase concentrations are history-dependent and differentially coupled to transport equations that contain spatial derivatives. Thus, in a system with n components and m phases, where every component can exist in every phase, LEA permits closure of the system with n partial differential equations, while a rigorous model for nonequilibrium requires mn. In groundwater modeling, historically n has been small, usually 3 or less, but increasing concerns with complex contaminants seem sure to stimulate interest in systems whose behavior demands many more components for accurate modeling. The petroleum industry has dealt successfully with systems of 10 or more components, but usually with IMPES techniques and always with LEA. It is highly unlikely that fully implicit treatment of mn equations, with n large and m around 3, will be practical in the foreseeable future. Hence, assuming that consideration of nonequilibrium is important, it will be critical to find simplifications that permit practical computations without unduly sacrificing accuracy.
If flux coefficients depend explicitly on phase concentrations, as in an IMPES
procedure, it should be relatively straightforward in a Newton-Raphson
iteration to set up and solve m flow
equations implicitly in a manner analogous to Young and Stephenson
[1983], using old values of phase concentrations, then solve 
transport
equations explicitly, including kinetics. This decouples flow from chemistry,
and caution may be called for in applying such an approach [ Valocchi and
Malmstead, 1992]. An intermediate possibility could be to somehow
approximate flowing phase concentrations (e.g., via local transient ordinary
differential equations, perhaps including data from adjacent cells),
short of full differential coupling,
given total component concentrations. The resulting n-equation model would
then resemble an equilibrium model, but with modified flux coefficients and
modified derivatives in the Newton-Raphson Jacobian. The approximate kinetic
techniques of Zaidel and Russo [1993] and Mackay et al. [1991],
or similar concepts, could be useful here. An adaptive implicit scheme
could enhance accuracy by allowing finer grids and time steps within practical
bounds; stability analysis would need to incorporate the kinetic relations.