This thesis analyzes partial differential equations related to the coupled surface and subsurface flows and develops efficient high order discontinuous Galerkin (DG) methods to solve them numerically. Specifically, the coupling of the Navier-Stokes and the Darcy's equations, which is encountered in the environmental problem of groundwater contamination through lakes and rivers, is considered. Predicting accurately the transport of contaminants by this coupled flow is of great importance for the remediation strategies. The first part of this thesis analyzes a weak formulation of the time-dependent Navier-Stokes equation coupled with the Darcy's equation through the Beavers-Joseph-Saffman condition. The analysis changes depending on whether the inertial forces are included in the interface conditions or not. The inclusion of the inertial forces (Model I) remedies the difficulty in the analysis caused by the nonlinear convection term; however, it does not reflect the physical interactions on the interface correctly. Hence, I also analyze the weak problem by omitting these forces (Model II) which complicates the analysis and necessitates an extra small data condition. For Model I, a fully discrete scheme based on the DG method and the Crank-Nicolson method is introduced. The convergence of the scheme is proven with optimal error estimates. The second part couples the surface flow and a convection-diffusion type transport with miscible displacement in the subsurface. Initially, I consider the coupled stationary Stokes and Darcy's equations for the flow and establish the existence of a weak solution. Next, imposing additional assumptions on the data, I extend the result to the nonlinear case where the surface flow is given by the Navier-Stokes equation. The analysis also applies to the particular case where the flow is loosely coupled to the transport, that is, the velocity field obtained from the flow is an input for the transport equation. The flow is discretized by combinations of the continuous finite element method and the DG method whereas the discretization of the transport is done by a combined DG and backward Euler methods. The scheme yields optimal error estimates and its robustness for fractured porous media is shown by a numerical example.