In the atmospheric surface layer (ASL), relationships between turbulent fluxes of momentum, heat and other scalars, and their mean gradients are described by Monin-Obukhov similarity theory (MOST), which states that any quantity of interest normalized by the proper scales (the height z, the kinematic surface shear stress u_*^2 = -tau_0/rho, the kinematic surface heat flux wt_0, and the buoyancy parameter g/Theta_0) will be a function only of the MOST stability parameter z/L. Although MOST is the standard framework in micrometerology, hydrology, ecology, and other disciplines used to interpret experimental measurements, to estimate turbulent fluxes (e.g. of trace gases) that cannot be calculated directly, and to parameterize turbulent transport in numerical weather and climate prediction models, a number of fundamental questions regarding the relationships between turbulent fluxes and their mean gradients in the ASL remain unanswered. When turbulent fluxes cannot be calculated directly, they are often estimated from MOST functions that relate the fluxes to mean gradients as a function of z/L. However, these functions are not predicted by theory and must be determined empirically. As a result, the connections between the observed behavior of these empirical curves and fundamental properties of turbulence in the ASL (e.g. the energy spectrum, turbulent kinetic energy budget, and integral lengthscales) are not well understood. The connections between the observed behavior of the dimensionless mean wind shear phi_m with z/L and observed properties of turbulence including the stability variation of the horizontal and vertical integral lengthscales are investigated using a theoretical framework coupled with ASL data from the Advection Horizontal Array Turbulence Study (AHATS). It is found that coefficients in empirical curves for phi_m(z/L) can be explained by accounting for the stability-dependence of the horizontal and vertical integral lengthscales, and the effect of buoyancy on the integral lengthscales is investigated in detail. When MOST functions are calculated from experimental data, a significant amount of scatter is typically found relative to empirical curves. The cause of this scatter is debated in the literature and, although many researchers attribute the scatter to random errors (i.e. due to an insufficient averaging period for the time average to converge to the true ensemble mean), others suggest that the scatter is physically meaningful and is due to physical processes that are neglected by MOST. In order to determine the extent to which random errors and physical processes each contribute to the observed scatter in experimental measurements, it is necessary to have a reliable method for estimating random errors. Most existing error estimation methods require an estimate of the integral time scale, which is highly sensitive to the method used to estimate it. A new method for estimating random errors, based on filtering the timeseries of the quantity of interest at multiple scales, is proposed and intercompared with existing error estimation methods. The new error estimation method is found to compare well with existing methods and does not require an a priori estimate of the integral timescale. In addition, it can be applied to estimate random errors in moments of any order. Using the new error estimation method, together with error propagation analysis, the influence of random errors on the observed scatter in MOST functions calculated from experimental data is explored. Errors in the turbulent fluxes and gradients are propagated to z/L, phi_m, and the dimensionless mean temperature gradient phi_h. Errors in the MOST stability parameter z/L are found to be large for unstable conditions, reaching values of 40% or more. Statistical hypothesis tests are used to demonstrate that random errors in z/L and phi_m or phi_h do not explain the observed scatter. In addition, deviations of phi_m from empirical curves are found to have strong diurnal variation and to increase with height, which suggests that physical processes related to the daytime growth of the convective boundary layer are responsible for the observed deviations from MOST. Finally, the large scale structure of convective organization in the unstable atmospheric boundary layer is investigated using large eddy simulation. Both observational and numerical studies have revealed that convective updrafts are organized into longitudinal rolls for slightly convective conditions and into open cells for highly convective conditions. However, the physical processes responsible for the transition from roll- to cellular- type convection, the nature of the transition (i.e. whether it occurs rapidly or gradually), and the criteria for the transition are not well understood. A suite of large eddy simulations spanning a range of -z_i/L, where z_i is the ABL depth and L is the Obukhov length, are conducted. Vertical profiles of the turbulent kinetic energy (TKE) anisotropy are found to collapse into three discrete regimes, which correspond to rolls, transitional structures, and cells in the velocity and temperature fields. The transitional regime is observed for 11