This study developed an ensemble four-dimensional variational (En4DVar) hybrid data assimilation (DA) system. Different from most of the available En4DVar systems that adopted ensemble Kalman Filter class or ensemble DA approaches to produce ensemble covariances for their hybrid background error covariances (BECs), it used a four-dimensional ensemble-variational (4DEnVar) system to obtain the ensemble covariance. The localization scheme for 4DEnVar applied orthogonal functions to decompose the correlation matrix so that it was implemented easily and rapidly. In terms of analysis quality and forecast skill, the En4DVar system was evaluated in the single-point observation experiments and observing system simulation experiments (OSSEs) with sounding and cloud-derived wind observations, using its standalone four-dimensional variational (4DVar) and 4DEnVar components as references. The single-point observation experiments visually verified the explicit flow-dependent characteristic of the BEC due to the introduction of the ensemble covariance from the 4DEnVar system. The OSSE-based sensitivity experiments revealed different contributions of the weight for the ensemble covariance in the En4DVar system to the forecasts in the Northern and Southern Extratropics and Tropics. A much higher weight for the ensemble covariance in a properly inflated hybrid covariance helped En4DVar produce the most reasonable analysis. The forecast initialized by En4DVar is overall better than by 4DVar and 4DEnVar, although the quality of En4DVar analysis is between those of 4DVar and 4DEnVar ensemble mean analyses. It indicates that the flow-dependent ensemble covariance provided by 4DEnVar dominantly contributes to the improvements in the En4DVar-initialized forecast, with certain but necessary constraint from the balanced climatological covariance.