Sobolev orthogonal polynomials of d variables on the product domain Omega := [a(1), b(1)] x ... x[a(d), b(d)] with respect to the inner product < f, g >(S) = (c)integral(Omega)del(k)f(x) . del(k)g(x)W(x)dx + Sigma(k-1)(i=m) lambda i del(i)g(p), k is an element of N, are constructed, where del(i) f, i = 0, 1, 2,..., kappa, is a column vector which contains all the partial derivatives of order i of f, x := (x(1), x(2), ..., x(d)) is an element of R-d, dx := dx(1)dx(2) ... dx(d), W(x) := w(1)(x(1))w(2)(x(2)) ... w(d)(x(d)) is a product weight function on O, wi is a weight function on [a(i), b(i)], i = 1, 2,..., d, lambda(i) > 0 for i = 0, 1,...,. k - 1, p = (p(1), p(2),..., p(d)) is a point in R-d, typically on the boundary of Omega, and c is the normalization constant of W. The main result consists of a generalization to several variables and higher order derivatives of some results which are presented in the literature of Sobolev orthogonal polynomials in two variables; namely, properties involving the integral part in <., .>(S), a connection formula, and a recursive relation for constructing iteratively the polynomials. To illustrate the main ideas, we present a new example for the HermiteHermite-Laguerre product weight function.