The abstract definition of the curvature tensor is good as is, see, if $\partial_\mu$ is the $\mu$th coordinate basis vector for some local chart, we have $\nabla_{\partial_\mu}\partial_\nu=\Gamma^\sigma_{\mu\nu}\partial_\sigma$, so we have $$ R(\partial_\mu,\partial_\nu)\partial_\rho=\nabla_{\partial_\mu}\nabla_{\partial_\nu}\partial_\rho-\nabla_{\partial_\nu}\nabla_{\partial_\mu}\partial_\rho=\nabla_{\partial_\mu}(\Gamma^\sigma_{\nu\rho}\partial_\sigma)-\nabla_{\partial_\nu}(\Gamma^\sigma_{\mu\rho}\partial_\sigma)=\partial_\mu\Gamma^\sigma_{\nu\rho}\partial_\sigma+\Gamma^\sigma_{\nu\rho}\nabla_{\partial_\mu}\partial_\sigma-(\mu\leftrightarrow\nu)=\\=\partial_\mu\Gamma^\sigma_{\nu\rho}\partial_\sigma+\Gamma^\sigma_{\nu\rho}\Gamma^\lambda_{\mu\sigma}\partial_\lambda-(\mu\leftrightarrow\nu)=(\partial_\mu\Gamma^\sigma_{\nu\rho}+\Gamma^\sigma_{\mu\lambda}\Gamma^\lambda_{\nu\rho}-(\mu\leftrightarrow\nu))\partial_\sigma, $$ so the formula $$ R(X,Y)Z=\nabla_X\nabla_YZ-\nabla_Y\nabla_XZ-\nabla_{[X,Y]}Z $$ gives the correct curvature formula in the case of torsion as well.
What confuses you is that in index notation, the expression $[\nabla_\mu,\nabla_\nu]$ actually corresponds to an antisymmetric expression made from second covariant derivatives, in abstract language, the torsionless curvature tensor can be expressed as $$ R(X,Y)Z=\nabla^2_{X,Y}Z-\nabla^2_{Y,X}Z, $$ and its THIS formula that needs to be modified for torsion. It's because $\nabla^2_{X,Y}Z$ corresponds to $X^\mu Y^\nu\nabla_\mu\nabla_\nu Z^\sigma$ in index notation. In abstract notation, $\nabla_X\nabla_Y Z\neq\nabla^2_{X,Y}Z$, because in the first term, $\nabla_X$ also acts on $Y$, while for the "second covariant derivatives", $Y$ is applied "from outside the differential operator".
In index notation, $$ \nabla_X\nabla_Y Z\Longleftrightarrow X^\mu\nabla_\mu(Y^\nu\nabla_\nu Z^\sigma) \\ \nabla^2_{X,Y}Z\Longleftrightarrow X^\mu Y^\nu\nabla_\mu\nabla_\nu Z^\sigma,$$ work out for yourself in component notation that these two are not equivalent.
The difference between the antisymmetric form of $\nabla_X\nabla_Y Z$ and $\nabla^2_{X,Y}Z$ will produce a term proportional to $\nabla_{\nabla_YX-\nabla_XY}$, which is $\nabla_{[X,Y]}$ iff $\nabla$ is torsionless, hence the need to modify $R(X,Y)Z=\nabla^2_{X,Y}Z-\nabla^2_{Y,X}Z$ if we have nonvanishing torsion, since the correct curvature tensor is given by $R(X,Y)Z=[\nabla_X,\nabla_Y]Z-\nabla_{[X,Y]}Z$.
Best Answer
Permit me the use of Latin indices instead of Greek indices and the convention $\nabla_a K_b=K_{b;a} $. So we wish to prove $\newcommand{\Tud}[3]{{#1}^{#2}_{\phantom{#2}{#3}}}$ $$\Tud{K}{a}{;b c} = \Tud{R}{a}{b c d} K^d$$ where $$\Tud{V}{a}{;b c} - \Tud{V}{a}{;c b} = \Tud{R}{a}{d c b} V^d$$ and $$K_{a ; b} + K_{b ; a} = 0$$
Differentiating the last equation, we get $$K_{a ; b c} + K_{b ; a c} = 0$$ so, relabelling and summing, $$K_{a ; b c} + K_{b ; a c} - K_{b ; c a} - K_{c ; b a} + K_{c ; a b} + K_{a ; c b} = 0$$ hence, $$K_{a; b c} + K_{a ; c b} = R_{b d a c} K^d + R_{c d a b} K^d$$ By the interchange symmetry $R_{a b c d} = R_{c d a b}$, and raising indices, we get $$\Tud{K}{a}{;b c} - \Tud{R}{a}{b c d} K^d = -(\Tud{K}{a}{;c b} - \Tud{R}{a}{c b d} K^d)$$
On the other hand, by the first Bianchi identity and antisymmetry, we have $$\Tud{R}{a}{d c b} = \Tud{R}{a}{b c d} + \Tud{R}{a}{c d b}$$ Hence we get $$\Tud{K}{a}{;b c} = \Tud{K}{a}{; c b} + \Tud{R}{a}{d c b} K^d = \Tud{K}{a}{;c b} + \Tud{R}{a}{b c d} K^d + \Tud{R}{a}{c d b} K^d$$ and therefore $$\Tud{K}{a}{;b c} - \Tud{R}{a}{b c d} K^d = \Tud{K}{a}{;c b} - \Tud{R}{a}{c b d} K^d$$ The conclusion follows.