It appears that you have some confusion regarding the basic notions of density operators and "dimension". Why $d^2$ dimensions "are required to describe" a density matrix isn't the right question to ask; density matrices are $d^2$ dimensional objects in the same sense that vectors in $\Bbb R^3$ are 3-dimensional (i.e., the cardinality of any basis set of $\Bbb R^3$ is 3). By convention, the dimensionality of a vector refers to the space of which the vector is a member. Or did you mean the minimum number of parameters needed specify a density matrix uniquely? More on that later.
One of the axioms of quantum mechanics is that to every quantum system there is associated a separable complex Hilbert space $(\mathcal{H}, +, \cdot, \langle \cdot|\cdot\rangle)$. The states of the system are all positive trace-class linear maps $\rho:\mathcal{H}\to \mathcal{H}$ for which $\mathrm{Tr}(\rho) = 1$. A quantum state$^\dagger$ $\rho$ is called a pure state if
$$\exists \psi\in \mathcal{H}: \forall \alpha \in \mathcal{H}: \rho(\alpha) = \frac{\langle \psi|\alpha \rangle}{\langle \psi|\psi\rangle}\psi.$$
To each pure state $\rho$ you can assign an element $\psi \in \mathcal{H}$ but the correspondence isn't one-to-one even if you restrict to pure states and impose the normalization condition (for instance, elements of the Hilbert varying by a global phase factor have essentially the same pure state). The take-home point is that it's wrong to refer to $\psi$ as the pure state itself; a correction perhaps would be that a pure state can be given by an element $\psi$ of the Hilbert space, only up to arbitrary rescaling.
Note that $\rho(\cdot)_\psi : \mathcal{H}\to\mathcal{H}$ is in fact the following linear operator
$$\rho(\cdot)_\psi = \frac{\langle \psi|\cdot\rangle}{\langle \psi|\psi\rangle}\psi = \frac{|\psi\rangle\langle\psi|}{\langle \psi|\psi\rangle}$$
which follows from the fact that $\mathcal H\otimes \mathcal H\simeq \mathcal H\otimes \mathcal H^*$ (cf. this answer). Once you choose an ordered basis, it can written down as a matrix as well.
Now this $\rho$ is what introductory textbooks call the "density operator" for a pure state, which is an endomorphism on $\mathcal{H}$ (a linear map from $\mathcal{H}$ to $\mathcal{H}$ and not merely an element of $\mathcal{H}$; this I believe was the root of your confusion). As you should know from your linear algebra classes: the dimension of space of linear maps between two vector spaces (say $V$ and $W$) is $\dim(V)\times \dim(W)$ (cf. proof), it is clear that the dimension of $\rho$ (w.r.t the space of such linear maps $\mathcal{H}\to\mathcal{H}$) is $\dim(H)\times \dim(H) = d \times d = d^2$. In infinite-dimensional cases things get a bit murkier but we leave that to the field of functional analysis.
Coming to mixed states, they're simply classical ensembles of pure states (encoding the lack of knowledge of the experimenter) and it can be shown from the spectral theorem that every mixed state can be represented as a non-unique countable convex combination of pure states $^{\dagger\dagger}$ i.e.,
$$\rho_{\text{mixed}} = \sum_i p_i\rho_{\psi_i}.$$
Introductory physics and quantum computing texts tend to write this as
$$\rho_{\text{mixed}} = \sum_i p_i|\psi_i\rangle\langle\psi_i|$$
with suitable normalization, but I do not find that notation particularly elegant.
Nevertheless, to answer your main question, linear operators $\mathcal{H}\to\mathcal{H}$ form a vector space structure $\mathcal{L}(\mathcal{H})$ and there exists a natural addition and scalar multiplication on that space. Clearly, the dimension is preserved upon a linear combination of $\dim^2(\mathcal{H})$ dimensional linear operators $\rho_{\psi_i}$ by the closure property of vector spaces (vector spaces are are closed under addition and scalar multiplication).
We need to be careful about the terminology here: when you say that the dimension of $\rho$ is $d^2$, it's actually referring to the dimension of $\mathcal{L}(\mathcal{H})$ i.e. the vector space of all linear maps from $\mathcal{\mathcal{H}}$ to $\mathcal{H}$ without any additional constraints. However, when you only allow convex combinations, the set of such $\rho$ s no longer has a vector space structure (cf. Moretti's answer); it's merely a convex set embedded in a vector space.
Furthermore, when you impose the conditions of positivity (self-adjointness, which is a basis-independent concept) and unit trace, the dimension (number of complex parameters) needed to specify $\rho$ then reduces to $$\underbrace{\frac{d^2-d}{2}}_{\text{off-diagonal terms}} + \underbrace{(d-1)}_{\text{diagonal terms}}.$$ That is, $(d^2-d)/2$ complex numbers are needed to specify the off-diagonal terms (the elements of a self-adjoint matrix are anti-symmetric about the diagonal) and only $(d-1)$ complex numbers are needed to specify the real diagonal terms (as the trace is constant). In this sense, the dimension of the "space" of density operators is $d(d+1)/2-1$ and not $d^2$. In colloquial language, it still makes sense to call it $d^2$ after all, because we usually call any $m\times n$ matrix as $mn$-dimensional, and you'll represent density operators as $d\times d$ matrices anyway, but it's useful to know the subtleties.
$^\dagger$: Not every quantum state can be represented by a density operator (cf. Chiral Anomaly's answer on Physics SE).
$^{\dagger\dagger}$: From Krein-Milman and Carathéodory's theorems it can be further shown that any state is a convex combination of at most $(\dim \mathcal H)^2$ pure states. Spectral theorem gives an even stronger result, namely that any state is a convex combination of at most $\dim \mathcal H$ pure states $|\psi_i\rangle\langle \psi_i|$ where $\psi_i$ are pairwise orthogonal unit vectors. For an elaborate discussion check the recent monograph titled Alice and Bob Meet Banach: The Interface of Asymptotic Geometric Analysis and Quantum Information Theory by Guillaume Aubrun and Stanisław J. Szarek (chapters 1-2).
