In the last few decades, laboratory experiments and direct numerical simulations of turbulent boundary layers, performed at low to moderate Reynolds numbers, have found very-large-scale motions (VLSMs) in the logarithmic and outer regions. The size of VLSMs was found to be 10-20 times as large as the boundary-layer thickness. Recently, few studies based on field experiments examined the presence of VLSMs in neutral atmospheric boundary-layer flows, which are invariably at very high Reynolds numbers. Very large scale structures similar to those observed in laboratory-scale experiments have been found and characterized. However, it is known that field measurements are more challenging than laboratory-based measurements, and can lack resolution and statistical convergence. Such challenges have implications on the robustness of the analysis, which may be further adversely affected by the use of Taylor’s hypothesis to convert time series to spatial data. We use large-eddy simulation (LES) to investigate VLSMs in atmospheric boundary-layer flows. In order to make sure that the largest flow structures are properly resolved, the horizontal domain size is chosen to be much larger than the standard domain size. It is shown that the contributions to the resolved turbulent kinetic energy and shear stress from VLSMs are significant. Therefore, the large computational domain adopted here is essential for the purpose of investigating VLSMs. The spatially coherent structures associated with VLSMs are characterized through flow visualization and statistical analysis. The instantaneous velocity fields in horizontal planes give evidence of streamwise-elongated flow structures of low-speed fluid with negative fluctuation of the streamwise velocity component, and which are flanked on either side by similarly elongated high-speed structures. The pre-multiplied power spectra and two-point correlations indicate that the scales of these streak-like structures are very large. These features are similar to those found in the logarithmic and outer regions of laboratory-scale boundary layers by direct numerical simulation and experiments conducted at low to moderate Reynolds numbers. The three-dimensional correlation map and conditional average of the three components of velocity further indicate that the low-speed and high-speed regions possess the same elongated ellipsoid-like structure, which is inclined upward along the streamwise direction, and they are accompanied by counter-rotating roll modes in the cross-section perpendicular to the streamwise direction. These results are in agreement with recent observations in the atmospheric surface layer. Furthermore, the influences of the Coriolis force, the domain size, and the subgrid-scale model on VLSMs are investigated.