A quantitative understanding of ground water flow characteristics in unconfined aquifers is important because of the prevalence of abstraction from, and pollution of these systems. The current understanding of ground water flow in unconfined aquifers is limited because of the dominance of horizontal flow modelling strategies used to represent unconfined flow processes. The application of horizontal flow principles leads to an ignorance of seepage-face formation and can not predict the complicated three-dimensional nature of the ground water flow that dominates at the ground water-surface water interface. This study aims to address some of these deficiencies by exploring the true three-dimensional nature of ground water flow including the formation of seepage faces at the ground water-surface water interface using numerical and laboratory techniques. A finite element model for simulating two-dimensional (vertical) variably saturated flow is developed and benchmarked against standard laboratory and field-scale solutions. The numerical features of the finite element model are explored and compared to a simple finite difference formulation. The comparison demonstrates how finite element formulations lead to a broader spatial averaging of material properties and a different method for the representation of specified flux boundaries. A detailed comparison analysis indicates that these differences in the finite element solution lead to an improved approximation to the partial differential equation governing two-dimensional (vertical) variably saturated flow. A laboratory analysis of unconfined ground water flow and associated solute transport characteristics was performed. The analysis focused upon unconfined flow towards a pumping well. The laboratory observations were reliably reproduced using a three-dimensional (axi-symmetric), variably saturated ground water flow model. The model was benchmarked against the ground water flow characteristics such as the seepage-face height and total flow rate. In addition, the model was shown to reliably reproduce the solute transport features such as travel times and streamline distributions. This is the first time that a numerical model has been used to reliably reproduce the solute transport characteristics near a seepage-face boundary where the three-dimensional flow effects are prevalent. The ability to reliably predict solute transport patterns in the seepage-face zone is important since this region is known to support vital microbially facilitated reactions that control nutrient cycling and contaminant attenuation. The three-dimensional travel time distribution near the seepage-face was compared to that predicted using a horizontal flow modelling approach derived from the basic Dupuit-Forchheimer equations. The Dupuit-Forchheimer based model indicated that horizontal flow modelling would under-estimate the total residence time near a seepage-face boundary, thereby introducing a considerable source of error in a solute transport analysis. For this analysis, a new analytical solution for the steady travel time distribution in an unconfined aquifer subject to a single pumping well was derived. The analytical model has identified, for the first time in the hydrogeology literature, the use of the imaginary error function. The imaginary error function is a standard transcendental function and an infinite series approach to evaluate the function was successfully proposed. The two-dimensional (vertical) ground water flow model was extended to handle the case where the flow is driven by density gradients near the ground water-surface water interface. The unsteady, two-dimensional, Galerkin finite element model of density-dependent ground water flow in variably saturated porous media is rigorously presented and partially benchmarked under fully saturated (confined) conditions. The partial benchmarking involved reproducing solutions to the standard Henry salt-water intrusion and the Elder salt-convection problems. The model was used in a standard density-coupled and a new density-uncoupled mode to elucidate the worthiness of the Henry and Elder problems as benchmark standards. A comparison of the coupled and uncoupled solutions indicates that the Henry salt-water intrusion problem has limited worthiness as a benchmark as the patterns of ground water flow are relatively insensitive to density-coupled effects. Alternatively, the Elder problem is completely dependent upon a correct representation of the density-coupled flow and solute transport processes. The coupled versus uncoupled comparison is proposed as a new test of the worthiness of benchmark standards. The Henry salt-water intrusion problem was further analysed in an attempt to alleviate some of the difficulties associated with this benchmark problem. The numerical model was tested against a re-evaluated version of Henry's semi-analytical solution for the coupled solute concentration distribution. The numerical model was used to propose a modified version of the Henry problem where the importance of density-coupled processes was increased. The modified problem was shown to have an improved worthiness as compared to the standard solution. The numerical model results were benchmarked against a new set of semi-analytical results for the modified problem. Certain advantages in using the modified problem as a test case for benchmarking the results of a numerical model of density-dependent ground water flow are identified. A numerical investigation of the patterns of density-driven ground water flow at the ground water-surface water interface was undertaken. The numerical model is shown to produce grid-independent results for a finely discretised domain. The pattern of discharge is controlled, in part, by two parameters. One describes the recharge applied to the aquifer, and the second describes the magnitude of the density differences between the fresh recharging fluid and the saline receiving fluid. The influence of dense intrusions upon the formation of seepage-face boundaries at the ground water-surface water interface under steady-state conditions was also investigated. Dense intrusions are shown to dominate the pattern of ground water flow only under mild recharge conditions, while seepage faces dominate the outflow pattern under strong recharge conditions. Therefore, the formation of seepage-face boundaries and dense intrusions are unlikely to coincide under the conditions examined in this study.