High-modulus asphalt mixtures are primarily employed to mitigate rutting distress in asphalt pavements. In addressing the effects of heterogeneous temperature fields under actual service conditions, a viscoelastic constitutive model for asphalt mixtures was derived, and a thermo-mechanically coupled viscoelastic dynamic analysis model incorporating temperature gradients was established. A quantitative method for characterizing depth-dependent viscoelastic properties was also proposed. Through dynamic modulus testing, broad-temperature-domain viscoelastic master curves were constructed using Boltzmann superposition principle and GaussAmp modeling. By integrating actual pavement temperature gradients, a unified 3D dynamic modulus surface model was developed. Results demonstrate that temperature gradients induce a significant stratified distribution of viscoelastic parameters: the dynamic modulus in the surface high-temperature zone (52 ℃) is reduced by 172% compared to that in deeper layers (40 ℃), and stress relaxation times are shortened by an order of magnitude. The developed pavement viscoelastic dynamic analysis model effectively characterizes the influences of temperature gradient, loading rate, and dwell time effects on viscoelastic responses of asphalt pavements. Thermo-coupled rutting analysis reveals that predicted rut depths under a non-uniform field deviate by 19.3% to 250.9% from those under a uniform temperature field, verifying that the effect of temperature gradient cannot be neglected. By utilizing a self-developed viscoelastic UMAT subroutine and a gradient parameter mapping method, a multiphysics-coupled analysis framework suitable for complex climatic conditions is established for high-modulus asphalt pavements, providing theoretical and technical support for rutting prediction and maintenance decision-making.