Abstract: Objectives: Ocean tide loading (OTL) is a significant error source affecting positioning accuracy for Global Navigation Satellite System (GNSS) stations. Currently, the systematic evaluation of ocean tide loading displacements (OTLD) derived from diverse OTL models on global-scale remain insufficient, while some International GNSS Service (IGS) analysis centers still only consider the 11 main ocean tide components, neglecting the impact of minor ocean tides (MOT). This gap hinders the optimization of OTL corrections for millimeter-level positioning and terrestrial reference frame (TRF) realization. The aim is to analyze threedimensional spatiotemporal characteristics of OTLD and inter-model discrepancies across global GNSS stations, with emphasis on incorporating MOT, using the six most recent representative OTL models. Methods: Six OTL models (FES2014b, TPXO9.5a, EOT20, GOT4.10c, NAO.99b, FES2004) were compared. Long-term, high-frequency OTLD time series (24 years, 15-minute sampling) at 122 globally distributed GNSS stations were generated using the HARDISP.F program, which considers 342 tidal constituents. Discrepancies were quantified via root mean square (RMS) and maximum (MAX) differences in vertical (U) and horizontal components, with analysis of spatiotemporal distributions and phase variations. Results: All models show strongest OTLD in coastal and island regions (peak vertical displacement of 85.74 mm). With exception of the Arctic Ocean, the Pacific, Atlantic, and Indian Oceans exhibit an overall trend of OTLD intensity attenuating as distance from shorelines increases, with distinct regional patterns: the Pacific Ocean region exhibits the most significant attenuation, with consistent OTLD time series features between its coastal and inland reference stations; the Atlantic follows with moderate attenuation, while the Indian Ocean shows the least pronounced attenuation; coastal stations in these two ocean regions display typical semi-diurnal tidal characteristics, whereas inland stations exhibit mixed tidal features, leading to notable differences in OTLD time series waveforms between coastal and inland stations. In contrast, the Arctic region demonstrates a distinctive non-monotonic attenuation pattern with lower decay rates, accompanied by significant inter-model discrepancies near coasts that converge at inland stations. Vertical OTLD (exceeding 40 mm at most coastal stations) is significantly larger than horizontal components (about 10 mm). Significant inter-model differences in magnitude and phase exist, particularly in high-latitude (e.g., Arctic coasts) and complex terrain regions (e.g., northeastern Pacific), among which the NAO.99b model shows a phase lag of approximately 6 hours for high tides. The RMS differences in vertical OTLD relative to FES2004 exceed 3 mm for most models, with a negative correlation between model differences and regional accuracy. NAO.99b shows the largest divergence (3.77 mm RMS vertically), with overall global accuracy far inferior to the other models. EOT20 and GOT4.10c follow the second, with bigger differences of up to 3.17 mm RMS in central Pacific islands and North Atlantic coasts. TPXO9.5a and FES2014b show minimal global differences, demonstrating relatively high accuracy worldwide overall, but with localized biases of up to 2.22 mm RMS vertically in regions like southeastern Indian Ocean and southwestern Atlantic. Conclusions: The results quantify critical OTLD spatiotemporal features and intermodel variabilities, demonstrating that one single OTL model would be insufficient for global high-precision GNSS applications, especially near coasts and islands. Our analysis underscores the importance of MOT and region-specific model performance. These findings would support developing optimized, regionally adaptive OTL correction strategies to enhance GNSS positioning accuracy and millimeter-level TRF realization.