For probability measures $μ,ν$ and $ρ$ define the cost functionals \begin{align*} C(μ,ρ):=\sup_{π\in Π(μ,ρ)} \int \langle x,y\rangle\, π(dx,dy),\quad C(ν,ρ):=\sup_{π\in Π(ν,ρ)} \int \langle x,y\rangle\, π(dx,dy), \end{align*} where $\langle\cdot, \cdot\rangle$ denotes the scalar product and $Π(\cdot,\cdot)$ is the set of couplings. We show that two probability measures $μ$ and $ν$ on $\mathbb{R}^d$ with finite first moments are in convex order (i.e. $μ\preceq_cν$) iff $C(μ,ρ)\le C(ν,ρ)$ holds for all probability measures $ρ$ on $\mathbb{R}^d$ with bounded support. This generalizes a result by Carlier. Our proof relies on a quantitative bound for the infimum of $\int f\,dν-\int f\,dμ$ over all $1$-Lipschitz functions $f$, which is obtained through optimal transport duality and Brenier's theorem. Building on this result, we derive new proofs of well-known one-dimensional characterizations of convex order. We also describe new computational methods for investigating convex order and applications to model-independent arbitrage strategies in mathematical finance.