The flow in the fully-developed wind-turbine array boundary layer is simulated using large-eddy simulation. The numerical framework is based on Fourier-spectral discretization in the horizontal and sixth-order compact finite differences in the vertical directions, with turbine forces modeled as actuator drag-disks. The numerical method is validated by reproducing previously published results of the flow in the atmospheric boundary layer (ABL) with and without wind turbines. The turbulent kinetic energy (TKE) budget in wind farms with increasing number of turbines in the streamwise direction is studied. The turbulent transport term is found to be negative above the rotor hub-height and positive below the hub-height, indicating a net transfer of TKE from above the hub-height to the lower half of the rotor region. Consistent with previous studies for an isolated turbine, the magnitudes of the turbulent transport, shear production and the TKE itself are found to increase with increasing number of turbines. An analytical model for the performance of infinitely large wind farms is evaluated using LES. A total of 15 simulations, covering aligned and staggered layouts, differing thrust coefficients, streamwise spacings of turbines and surface roughness heights of the ABL, are considered. The model is found to provide correct trends for the LES data, but is found to consistently underpredict the power and thrust coefficients, by approximately 5-20% and 1-12%, respectively. The under-prediction is partly explained by the assumptions of inviscid flow and validity of the classical actuator-disk theory, inherent in the one-dimensional model. The LES results suggest that the empirical parameter employed by the analytical model can be estimated based on the turbine loading.