A new multivariate Archimedean copula estimation method is proposed in a non-parametric setting. The method uses the so called Geometrically Designed splines (GeD splines), recently introduced by Kaishev et al. (2006 a,b) [10] and [11], to represent the cdf of a random variable W θ , obtained through the probability integral transform of an Archimedean copula with parameter θ. Sufficient conditions for the GeD spline estimator to posses the properties of the underlying theoretical cdf, K(θ, t), of W θ , are given. The latter conditions allow for defining a three-step estimation procedure for solving the resulting non-linear regression problem with linear inequality constraints. In the proposed procedure, finding the number and location of the knots and the coefficients of the unconstrained GeD spline estimator and solving the constraint least-squares optimisation problem, are separated. Thus, the resulting spline estimatorK(θ, t) is used to recover the generator and the related Archimedean copula by solving an ordinary differential equation. The proposed method is truly multivariate, it brings about numerical efficiency and as a result can be applied with large volumes of data and for dimensions d ≥ 2, as illustrated by the numerical examples presented.