Laminar separation bubbles on a flat plate boundary layer in the presence of free-stream turbulence (FST) were investigated by means of Direct Numerical Simulations (DNS). A suction/blowing velocity distribution was applied along the free-stream boundary of the computational domain to induce separation on the flat plate. For numerically generating free-stream turbulence, isotropic grid turbulence, which is obtained from a superposition of eigenmodes from the continuous spectrum of Orr- Sommerfeld and Squire Equations, was introduced at the inflow boundary. The effect of the spanwise extent of the computational domain was investigated by carrying out computations with two different spanwise domain widths. The main characteristics of the separation bubble, such as the bubble length and the skin-friction distribution were very similar for both spanwise domain sizes. However, for the wider domain the 2-D "rollers" were modulated in the spanwise direction and broke up earlier than for the narrow domain. For the narrow domain counter-rotating streamwise vortices appeared ("braids"). Detailed numerical simulations were performed to investigate the effect of the free-stream turbulence energy spectrum. It was found that the transition location was essentially independent of the integral length scale of the free-stream turbulence. Also, the dependence of the separation length on the integral length scale was found to be very weak for the range of length scales considered in our studies. In contrast, the separation length was significantly reduced already for relatively low free-stream turbulence intensity (0.1%) when compared to the baseline case with zero FST. When the FST was increased further the length and height of the bubble continued to decrease. Instantaneous flow field visualizations revealed that the spanwise coherence of the dominant 2D structures was weakened with increasing FST intensity. In addition, proper orthogonal decomposition (POD) analyses of the instantaneous flow data revealed that streamwise "vortical" structures became dominant for high FST levels. Based on a detailed analysis of the time-dependent flow field and a comparison between DNS results and linear stability theory (LST) calculations, it was found that for free-stream turbulence intensities up to 2%, transition in the bubble was still due to an inviscid (Kelvin-Helmholtz) instability of the inflectional velocity profile in the separated flow region, and not due to nonlinear bypass mechanisms.