We investigate, by numerical lattice simulations, the static quark-antiquark potential, the flux tube properties and the chiral condensate for $N_f = 2+1$ QCD with physical quark masses in the presence of strong magnetic fields, going up to $eB = 9$ GeV$^2$, with continuum extrapolated results. The string tension for quark-antiquark separations longitudinal to the magnetic field is suppressed by one order of magnitude at the largest explored magnetic field with respect to its value at zero magnetic background, but is still non-vanishing; in the transverse direction, instead, the string tension is enhanced but seems to reach a saturation at around 50 % of its value at $B = 0$. The flux tube shows a consistent suppression/enhancement of the overall amplitude, with mild modifications of its profile. Finally, we observe magnetic catalysis in the whole range of explored fields with a behavior compatible with a lowest Landau level approximation, in particular with a linear dependence of the chiral condensate on $B$ which is in agreement, within errors, with that already observed for $eB sim 1$ GeV$^2$.