Don't trust the source power spectrum results

re-generated BEM surfaces

Have you checked visually the surfaces?
Open the 3 BEM layers and the low-resolution cortex in the same figure, set a different color for the inner skull for improved readability post the screen capture here.

re-calculated the OpenMEEG BEM head model

  • Have you seen any warning in a message box or in the Matlab command window that indicates that some dipoles are not positioned correctly?
  • For debugging further the OpenMEEG solution: right-click on the forward model > View lead field vectors, navigate between sensors and see if there are any dipoles that do not follow the trend of their neighbors.
  • I suggest you try comparing the OpenMEEG BEM results with a spherical model: compare the source maps, and compare the leadfields (select both the BEM and spherical models, right-click > View lead field vectors)

So, it seems that the BEM head model is not right due to a problem with the Freesurfer surfaces. I then recalculated the surfaces using CAT12 and ran through the steps above. However, the results were very similar to using the Freesurfer surfaces, i.e. not good.

The quality of the BEM solution depends mostly on the BEM surfaces, which are generated by Brainstorm. Whether you use CAT12 or FreeSurfer to generate the cortex surface influences very little these BEM layers.
It is possible that one cortex surface (FreeSurfer or CAT) has more vertices outside the inner skull (computed by Brainstorm) than the other, which could cause instabilities. But this is something you can fix with a right-click on the cortex > Force inside skull.

computed sources using Current Density and sLORETA maps

The minimum norm maps without normalization (current density maps) look similar and OK for both the individual head model and the template head model. The problem is probably not in the surfaces, but rather in the sLORETA normalization (widely used in the EEG literature, but not much promoted by Brainstorm and MNE-Python developers). I suggested you try with dSPM instead, as done in this tutorial:
https://neuroimage.usc.edu/brainstorm/Tutorials/RestingOmega#Source_estimation

Are these recordings only one task (eg. resting or sleeping)? Otherwise, if you have two states (eg. a task vs. rest), compute the noise covariance from the resting periods, or normalize the source maps with a Z-score wrt the rest recordings.

Reference Issue. The EEG data was already referenced to the average. Re-referencing again to the average made no difference.

In the forum thread I pointed out, the problem was related with the type of some data channels (adding or removing the mastoid electrodes to the group of EEG channels), but this is most likely not a problem here, because the recordings were converted to avg ref before being imported in Brainstorm.

I could not compute noise covariance, when right-clicking on the channel file > Noise covariance, the options are: no noise modeling, or import from file or from matlab.

If you want to compute the noise covariance matrix from recordings, you need to right-click on the RECORDINGS, not the channel file.