This study develops a physics constrained Gaussian process (GP) constitutive model for rate dependent creep and demonstrates its integration within a nonlinear finite element framework. The proposed approach provides a nonparametric alternative to conventional Norton type power law models while preserving thermodynamic admissibility. The model is constructed using uniaxial constant load creep data for 9Cr–1Mo steel at 550℃ across seven stress levels. A hierarchical two stage learning strategy is adopted. First, creep strain–time responses at each stress level are analyzed independently using GP regression, where monotonicity is imposed through a probit likelihood [1], rendering the posterior analytically intractable. Expectation propagation [2] is thus employed to obtain an efficient Gaussian approximation. This strategy enables stable inference of creep strain rates without relying on numerically sensitive differentiation of experimental measurements. The constitutive representation is restricted to the primary and secondary regimes, excluding tertiary acceleration associated with damage evolution. In the second stage, a strain hardening rate relation is learned directly from the inferred rates. Accumulated inelastic strain is employed as a phenomenological surrogate for microstructural evolution, allowing the model to capture the transition from primary hardening toward steady state creep. Strict positivity of the creep rate is enforced through an exponential mapping, ensuring satisfaction of the positive dissipation requirement imposed by the second law of thermodynamics. Additional physical constraints are incorporated at virtual observation points using expectation propagation. Positive stress sensitivity is enforced across all regimes, while a non positive strain rate sensitivity is imposed only within the primary creep regime to reflect strain hardening. For finite element implementation, an implicit time integration algorithm with a consistent tangent modulus is developed. The GP model supplies both the creep rate and its derivatives with respect to stress and strain, enabling Newton iterations at the material point level. The proposed framework establishes a direct pathway from experimental data to finite element simulation, enabling flexible representation of rate dependent behavior while maintaining thermodynamic consistency and uncertainty quantification.