Why averaging in Spectral magnitude domain not in complex domain to estimate spectrum of a process
You seem to understand the issues at hand.  In your situation, if the trial runs would be identical, except for the noise, averaging the signals (or the complex values) will definitely work, and I recommend that it how you do it.

The notion of averaging magnitudes (or magnitudes squared) across DFTs comes from when you are chunking a signal into multiple DFTs and want to aggregate the results.  In that case, if you average the raw DFTs you will get all sorts of cancellation/reinforcement effects depending on the frequency values of the signal compared to your DFT frame length.  When you average the magnitudes, the phase values are ignored and you will get more desirable results.

As is often the case wiht policy, it sometimes gets followed blindly without the understanding of the rationale behind the policy.  That may be what is behind some of the recommendations you have heard.

If you are trying to decompose your signal (assuming it is periodic) into its constituent harmonics, retaining the phase information is essential.

Hope this helps,

Ced