<abstract><p>The integration of distributed energy resources (DERs) and, therefore, power electronic devices into distribution networks leads to harmonic distortion injection. However, studying harmonic distortion solely through deterministic approaches presents challenges due to the inherent random behavior of DERs. This study introduced a strategy that leverages PowerFactory's harmonic load flow tool. By combining it with Python co-simulation, probabilistic load flows can be developed. These load flows utilize current sources to represent harmonic distortion emitters with predefined harmonic spectra. The proposed strategy was implemented on a real network, where two different capacities of DERs were integrated at various locations within the network. The distributions for the total harmonic distortion of voltage ($ THD_{v} $) and the total harmonic distortion of current ($ THD_{i} $) were obtained 24 hours a day in nodes and lines of the network. The procedure allowed considering the uncertainty associated to the DERs integration in distribution networks in the study of harmonic distortion, which, speaking from a simulation approach, is scarce in the literature.</p></abstract>