Assessing　the　seismic　capacity　of　pier　columns　is　a　crucial　element　in　the　performance　-based　seismic　design　of　bridges.　Such　assessment　necessitates　a　probabilistic　approach　to　accurately　determine　the　marginal　probability　distributions　of　seismic　capacities　and　to　characterize　the　dependencies　among　these　variables.　In　response　to　this　need,　this　paper　employs　Multi　-Output　Gaussian　Process　Regression　(MO-GPR),　a　probabilistic　machine　learning　method,　to　jointly　model　the　multiple　seismic　capacities　of　pier　columns.　We　initially　introduce　a　probabilistic　seismic　capacity　model　that　utilizes　MO-GPR　for　pier　columns　and　validate　its　predictive　accuracy　in　comparison　to　Bayesian　linear　regression　and　existing　empirical　methods.　Subsequently,　the　methodology　is　augmented　by　the　integration　of　hierarchical　modeling　within　the　MO-GPR　framework,　resulting　in　a　Multi　-Output　Hierarchical　Gaussian　Process　Regression　(MO-HGPR)　model　that　effectively　estimates　intraclass　correlation　coefficients　for　specific　types　of　datasets.　It　is　postulated　that　these　correlation　coefficients　also　reflect　correlations　associated　with　multiple　components　of　the　real　structure.　This　study　employs　MO-HGPR　and　MO-GPR　separately　to　investigate　the　potential　correlations　of　seismic　capacities　among　pier　columns　and　distinct　limit　states.　The　results　demonstrate　that　the　MO-GPR　model　exhibits　superior　prediction　accuracy　and　more　effectively　portrays　uncertainty　compared　to　existing　empirical　models.　Moreover,　the　correlations　of　seismic　capacities　among　piers　and　limit　states　are　both　robust　and　significantly　impact　the　seismic　fragility　of　bridges.　This　finding　highlights　the　essential　nature　of　considering　capacities　correlations　in　seismic　fragility　or　risk　assessment　processes.