$\newcommand{\R}{\mathbb{R}}$
Welcome to MathOverflow! My conjecture is as follows. Let $a_1,\dots,a_N$ be the distinct points in question. For $i\in[N]:=\{1,\dots,N\}$ and each $k=(k_1,\dots,k_N)\in\R^N$, let
\begin{equation}
X_i(k):=\{x\in[0,1]^d\colon|x-a_i|^p-|x-a_j|^p\le k_i-k_j\ \ \forall j\in[N]\setminus\{i\}\}. \tag{1}
\end{equation}
Note that, if $k_i=0$ for all $i$, then the family $X(k):=(X_i(k))_{i\in[N]}$ is the Voronoi tesselation for the points $a_1,\dots,a_N$. So, one may refer to $X(k)$ in general as the $k$-Voronoi tesselation.
Conjecture 1 For some $k\in\R^N$, the cells $X_i(k)$ of the $k$-Voronoi tesselation $X(k)$ are all of the same $d$-volume, $1/N$.
Let us denote such a vector $k$ by $k_*$.
Conjecture 2 The optimal transportation of the uniform distribution on the set $\{a_1,\dots,a_N\}$ to the uniform distribution on the $d$-cube $[0,1]^d$ is given by the transportation of the $\frac1N$-mass at each point $a_i$ to $\frac1N\,\times\big(\text{the uniform distribution on the cell }X_i(k_*)\big)$.
So, the $p$th power of the $p$-Wasserstein distance will be
$$\sum_{i=1}^N\int_{X_i(k_*)}|x-a_i|^p\,dx.
$$
Informal justification: Let $m_i(A)$ denote the mass transported from a point $a_i$ to a Borel set $A\subseteq[0,1]^d$. We have to minimize
\begin{equation*}
\sum_i\int_{[0,1]^d}|x-a_i|^p m_i(dx)
\end{equation*}
given that $m_i\ge0$, $\int_{[0,1]^d}m_i(dx)=1$ for all $i$, and $\sum_i m_i(dx)=dx$. Varying the measures $m_i$ and using Lagrange multipliers, we have
$|x-a_i|^p=k_i+\mu(x)$ for some $k=(k_1,\dots,k_N)\in\R^N$, some function $\mu$, all $x$, and all $i$ such that $x$ is in the support set (say $S_i$) of the measure $m_i$. It follows that $|x-a_i|^p-|x-a_j|^p=k_i-k_j$ for all $x$ and all $i,j$ such that $x\in S_i\cap S_j$. This gives rise to formula (1).
