For probability measures mu, nu, and rho, define the cost functionals C(mu, rho) := sup(pi is an element of Pi(mu, rho)) integral < x, y >pi(dx, dy) and C(nu, rho) := sup(pi is an element of Pi(nu, rho)) integral < x, y >pi(dx, dy), where <center dot,center dot > denotes the scalar product and Pi(center dot,center dot) is the set of couplings. We show that two probability measures mu and nu on R-d with finite first moments are in convex order (i.e., mu <=(c) nu) iff C(mu, rho) <= C(nu, rho) holds for all probability measures rho on R-d with bounded support. This generalizes a result by Carlier. Our proof relies on a quantitative bound for the infimum of integral fd nu - integral fd mu over all 1-Lipschitz functions f, which is obtained through optimal transport (OT) duality and the characterization result of OT (couplings) by Ruschendorf, by Rachev, and by Brenier. 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.