In this paper, we investigate massive multi-input multi-output (MIMO) high frequency (HF) skywave communications. We first introduce a model for HF skywave massive MIMO channels within the orthogonal frequency division multiplexing transmission framework by using the matrix of sampled steering vectors. The steering vectors vary across different subcarriers due to the effect of the propagation delay across the largescale antenna array. Specifically, we derive a wideband beam based channel model and show that the beam domain statistical channel state information (CSI) is frequency-independent. Then, we consider minimum mean-squared error based uplink receiver and downlink precoder with perfect CSI at the base station (BS). With a large number of antennas at the BS, the sum-rate can be asymptotically increased proportionally to the number of user terminals (UTs) while the transmit power per UT is scaled down inverse-proportionally to the number of antennas. Simulation results demonstrate very significant performance advantages of the proposed HF skywave massive MIMO system.