Ok, correct me if I'm wrong: you are claiming that a time invariant, non linear, two port network, having a 3rd degree polynomial transfer function, excited with a pure sine at the input port, shows not only the 2nd and 3rd harmonics at the output, but the fundamental H1 will also change from what could be a strictly proportional ratio, to the input sine amplitude. From this behavior you conclude that THD, defined as the geometric average of harmonics amplitudes to the fundamental not entirely consistent.
Your claim is correct, and you could intuitively look at such a system as having IMD, and then the H3-H2 component goes to the H1 bin, which obviously always happens. In fact, if you extend the ugly trigonometric algebra, you will mathematically get this H3-H2 intermodulation product too. Ok, so what? This second order intermodulation (adding/subtracting to H1, depending on the phase) products will always be much smaller that the first order intermodulation products (H2, H3) and the linear fundamental H1. Therefore, except for pathological cases where the distortion levels are very high, these 2nd (and 3rd... etc... depending on the non linear transfer function) can usually be safely ignored. It would be a good exercise for you to work up a practical example of an HIFI audio system where these high order intermodulation products would significantly alter the calculated THD. I'm sure characterizing the famous $250K Wavac amplifier using THD is incorrect
THD is nothing but a definition, and the number provided is just that: a synthetic number. The limitations of this definition is the reason why (for example) DSL signals (where noise and distortions are critical, since they both impact the channel bandwidth and the adjacent channel SNR) are not characterized using THD, but other metrics are used. In fact, I am myself not aware of another domain, other than audio, where the THD metric is consistently used.