Why does not ContourPlot provide double roots?

ContourPlot[(x - 1)^2 (y - 3)^2,
 {x, 0, 5}, {y, 0, 6},
 Contours -> 10.0^Range[-2, 1],
 PlotPoints -> 50,
 Exclusions -> {x == 1, y == 3},
 ExclusionsStyle -> Directive[AbsoluteThickness[1/4], Red]]

An alternative approach:RegionPlot + ImplicitRegion

RegionPlot @ ImplicitRegion[(x - 1)^2 (y - 3)^2 == 0, {{x, 0, 5}, {y, 0, 6}}]