A numerical model is presented to describe solute transport in groundwater coupled to sorption by plant roots, translocation into plant stems, and finally evapotranspiration. The conceptual model takes into account both Root Concentration Factor, RCF, and Transpiration Stream Concentration Factor, TSCF for chemicals which are a function of Kow. A similar technique used to simulate the solute transport in groundwater to simulate sorption and plant uptake is used. The mathematical equation is solved using finite difference technique to solve for the concentration at any grid cell with respect to time. The new package is integrated into SEAM3D to create a new SEAM3D Plant Uptake Package, or PUP. The model is then verified by comparing results for root sorption in one side to the SEAM3D Reaction Package, and results for plant uptake to the SEAM3D Source Sink Mixing Package. The verification results showed an excellent match, which led to using the new package in a series of design application scenarios to evaluate phytoremediation effect. Hypothetical design scenarios included: 1) the effect of a phytoremediation system dimensions, 2) the effect of phytoremediation plant density or maximum ET rate, 3) the effect of out-flux of the phytoremediation with respect to the natural aquifer in-flux, and 4) the effect of using a phytoremediation system when the source of contamination is removed. For all the previous study cases, the results evaluate the effect on: 1) contaminant concentrations downstream the source (expressed in plume length at a concentration 1% of the source concentration), 2) solute mass removal from the aquifer, and 3) mass-flux changes at different cross-sections downstream the contaminant source. The results indicating the followings: 1) the width of the phytoremediation system, WET, has a limited effect on the solute mass-removal; 2) high tree density close to the contaminant source has a greater effect on solute mass removal relative to uniform density of trees planted over the entire plume; 3) the width of the ET area will have only a slight effect on the mass removal if the TSCF value is small; 4) as the value of TSCF gets lower, the efficiency of solute mass uptake is lower, and thus the solute concentration in groundwater is higher regardless of the quantity of water transpired; 5) dynamic steady-state plume dimensions (specially the plume length) are affected by the groundwater in-flux, which will control the dimensions and density of a phyto system; 6) splitting the phyto system into two halves does not have the same outcome of having one piece of area closer to the contamination site; 7) using a phyto system after the contamination source is removed led to increasing the solute concentration in the areas of the trees and decreases the concentration in the areas downstream the trees. The alternative model gives more options for simulation of solute mass uptake by plants by making use of field and lab data between the solute dissolved concentration in groundwater C, and solute mass in tree's core M to select a modeling category of three: Linear (ISO-1), Freundlich (ISO-2), and Langmuir (ISO-3). Each modeling option depends on the designer selection according to the fitted equation parameters between, C and, M. In terms of conservative results, ISO-1, and ISO-2 give less mass removal results than ISO-3 in case of sources with low concentrations. ISO-2, and ISO-3 give less mass removal results than ISO-1 in case of sources with high concentrations. / Ph. D.