The global energy transition that is shifting the industry to renewable energy sources is at the center of many countries' strategic plans as demonstrated by the increasing number of offshore wind farms connected to the grid. However, modern-day wind turbine with full-scale back-to-back converter is considered as a main source of harmonics, which can be amplified through the network due to low frequency resonance points created by the combination of long capacitive HVAC submarine cables and the many inductive transformers. The amplification of harmonics could lead to a grid non-compliance, and hence the developers of Offshore Wind Power Plants (OWPPs) might be obligated to install inflexible and expensive filters. Prior to the construction of a wind farm and during the design stage, it is necessary to perform harmonic studies to estimate the harmonic emission from the wind farm and to predict whether there would be a harmonic non-compliance problem and the necessity or not to plan for remedial measures. The industry currently utilizes the IEC second summation rule to estimate harmonic distortion at the PCC of OWPPs. However, there is many proof in the literature that conclude the ineffectiveness of the summation rule as it leads to inaccurate estimation of harmonics compared to actual measurements. Therefore, there is a rising necessity to develop new accurate methods to estimate the harmonic distortion at the PCC of OWPPs. It has been shown in the literature that harmonics injected by wind turbine converters are statistically random, and both of their magnitudes and phase angles can be represented by probability distribution functions. In contrast, the IEC summation rule utilizes a deterministic 95th percentile magnitude and indirectly compensates for the phase angle by an exponent. These two assumptions could be the main sources of its inaccuracy. Therefore, the objective of this thesis is to develop a new method to estimate harmonic distortion at the PCC of an OWPP using probability theory considering the probability distribution functions of harmonics. First, a literature review is performed to acquire probability distribution functions of different harmonics based on actual measurements of a full-scale back-to-back converter, and on the modelling of the converter and other elements of the wind farm for harmonic studies. Then, an offshore wind power plant is modelled in DIgSilent PowerFactory software to perform both phase-correct (proposed method) and IEC harmonic load flows. The harmonic load flow for the proposed method is automated using Monte Carlo programmed in Python codes to calculate the harmonic distortions at the PCC, from which histograms can be computed to represent the statistical nature of the results. As the 95th percentile of the harmonic distortion at the PCC is the value used in harmonic studies to prove compliance, it is computed and used for comparison to the IEC summation rule results for four case studies. In case A, the IEC second summation rule underestimates low order harmonics such as the 5th and 7th harmonics and overestimates the 8th harmonic compared to the Monte Carlo method. It is demonstrated that the probability distribution functions of the magnitude and phase angle of the harmonics significantly impact the distribution of the harmonic distortion at the PCC. It is shown in the thesis that the probabilistic phase-correct (Monte Carlo) method provides closer results to the expected results, as seen from measurements, than IEC does, hence with the Monte Carlo method the uncertainty in grid-compliance studies is lower. In case B, it is studied whether the variation of the grid impedance could lead to grid code violation that could not have been predicted using a fixed grid impedance. The variation of the grid impedance reveals a big difference between the lowest and highest harmonic distortion that could be as high as 5 times, while the error between the probabilistic and IEC summation rule methods could be as twice high. This case concludes that the variation of the grid impedance is an important factor that should not be neglected in harmonic studies. In case C, the nonlinear behavior of the converter is included in the model and the impact of the frequency coupling phenomenon on the harmonic distortion at the PCC is studied. The phenomenon results in low errors when compared to the probabilistic method, however the error is quite high when compared to IEC summation rule. In case D, the frequency coupling phenomenon is analyzed while varying the grid impedance, which results in higher errors when compared to IEC summation rule. The results of the last two cases support the argument to include the frequency coupling phenomenon in harmonic studies. Therefore, it is concluded that the accuracy of harmonic studies could be improved by including the statistical nature of the wind turbine harmonics, the variation of the grid impedance and the frequency coupling phenomenon of the wind turbine converters.