I think your (7) is wrong, perhaps you didn't apply the product rule? Anyway, here's what I got (I'm going to call the dimension $n$ rather than $d$ because I find it too confusing with the $d$ from the differentials lying around):
\begin{align}
g(r)&=r^n\int_{B_1(0)}f(x+rz)\,dz
\end{align}
So,
\begin{align}
g'(r)&=nr^{n-1}\int_{B_1(0)}f(x+rz)\,dz + r^n\int_{B_1(0)}\langle\nabla f(x+rz), z\rangle\,dz\\
&=\frac{n}{r}\int_{B_r(x)}f(y)\,dy + \int_{B_r(x)}\left\langle \nabla f(y),\frac{y-x}{r}\right\rangle\,dy\\
&=\frac{n}{r}\int_{B_r(x)}f(y)\,dy + \frac{1}{r}\int_{B_r(x)}\sum_{i=1}^n\frac{\partial f}{\partial y^i}(y)\cdot (y^i-x^i)\,dy\\
&=\frac{n}{r}\int_{B_r(x)}f(y)\,dy+\frac{1}{r}\int_{B_r(x)}\sum_{i=1}^n\bigg[\frac{\partial }{\partial y^i}\left(f(y)\cdot (y^i-x^i)\right)-f(y)\bigg]\,dy
\end{align}
Note that the first integral and the last term from the summation cancel out. What remains is the divergence of the vector field $y\mapsto f(y)\cdot (y-x)$, so after applying the divergence theorem, we get
\begin{align}
g'(r)&=\frac{1}{r}\int_{\partial B_r(x)}\sum_{i=1}^n[f(y)\cdot (y^i-x^i)]\cdot \nu^i\,d\sigma_r(y)\\
&=\int_{\partial B_r(x)}f(y)\,d\sigma_r(y),
\end{align}
where the last equality is because $\frac{y-x}{r}=\nu$ is the unit outward normal, so the inner product appearing above is just $1$. This is how to do the calculation if we assume say $f$ is $C^1$ (so that all the steps above are justified).
If you do not wish to assume so much regularity on $f$, then we can apply the usual tricks of analysis: replace $f$ by a smoothed out version by convolving against a $C^{\infty}_c$ function, apply the previous result, then take limits to deduce the original case. More explicitly, let $\phi\in C^{\infty}_c(\Bbb{R}^n)$ be such that $\int_{\Bbb{R}^n}\phi\,d\lambda=1$, and for each $\epsilon>0$, define $\phi_{\epsilon}(y)=\frac{1}{\epsilon^n}\phi\left(\frac{y}{\epsilon}\right)$.
Since the differentiability of $g$ is a pointwise matter, we may modify $f$ outside some large ball so that $f$ is continuous and has compact support, and this will not affect the value of the integral defining $g$. Then, by the standard arguments regarding convolutions (which you can find in any decent textbook, eg Folland), $f_{\epsilon}:=f*\phi_{\epsilon}$ will be smooth, and as $\epsilon\to 0^+$, we have $f_{\epsilon}\to f$ uniformly on compact subsets of $\Bbb{R}^n$. Now, define
\begin{align}
g_{\epsilon}(r):=\int_{B_r(x)} f_{\epsilon}\,d\lambda
\end{align}
Then, by the previously established case, we know $g_{\epsilon}$ is differentiable and that
\begin{align}
g_{\epsilon}'(r)&=\int_{\partial B_r(x)}f_{\epsilon}\,d\sigma_r
\end{align}
Since $f_{\epsilon}\to f$ uniformly on compact subsets of $\Bbb{R}^n$ it follows that $g_{\epsilon}\to g$ uniformly on compact subsets of $\Bbb{R}$ and that $g_{\epsilon}'(\cdot)\to \int_{\partial B_{(\cdot)}(x)}f\,d\sigma_{(\cdot)}$ uniformly on compact subsets of $\Bbb{R}$. As a result of this uniform convergence, we can justify interchanging the limit $\epsilon\to 0^+$ with the derivative to deduce that $g$ is differentiable and that has the derivative we expect
\begin{align}
g'(r)&=\int_{\partial B_r(x)}f\,d\sigma_r.
\end{align}
The other approach to solving this question is to first establish the identity
\begin{align}
g(r)&:=\int_{B_r(x)}f(y)\,dy\\
&=\int_0^r\int_{\partial B_t(x)}f(y)\,d\sigma_t(y)\,dt,
\end{align}
From here, continuity of $f$ will allow us to use the standard single-variable fundamental theorem of calculus to deduce the answer. In this case, I think using spherical coordinates parametrization is the quickest way of getting the answer; but of course one should be comfortable with manipulating determinants, and the definition of integration on manifolds using charts, and changing variables.
Let $\xi:U\subset \Bbb{R}^{n-1}\to S^{n-1}=\partial B_1(0)\subset\Bbb{R}^n$ be a spherical coordinate parametrization. By this I mean a mapping $(\theta^1,\dots, \theta^{n-1})\mapsto \xi(\theta^1,\dots, \theta^{n-1})$ of an open subset $U$ of $\Bbb{R}^{n-1}$ onto an open subset of the unit sphere $S^{n-1}$, which covers "most" of it in the sense that the measure of the difference $\sigma_1\left[S^{n-1}\setminus \xi(U)\right]=0$. In two dimensions this is just $\theta\mapsto (\cos\theta,\sin\theta)$, in $3$ dimensions there are two angles, and so on. In higher dimensions, there is certainly an explicit formula, but that's not relevant. Let $\zeta:(0,\infty)\times U\to \Bbb{R}^n\setminus \{x\}$ be the mapping $(r,\theta)\mapsto x+r\cdot\xi(\theta)$. Then,
\begin{align}
g(r)&:=\int_{B_r(x)}f(y)\,dy\\
&=\int_{(0,r)\times U}f(\zeta(t,\theta))\cdot\left|\det D\zeta_{(t,\theta)}\right|\, d(t,\theta)\\
&=\int_0^r\int_Uf(x+t\xi(\theta))\cdot \sqrt{\left|\det \bigg[\bigg\langle D\zeta_{(t,\theta)}(e_i), D\zeta_{(t,\theta)}(e_j)\bigg\rangle\bigg]_{i,j\in\{1,\dots, n\}}\right|}\,d\theta\,dt.
\end{align}
Now, we should simplify that determinant. Observe first that the matrix representation of $D\zeta_{(t,\theta)}$ (relative to the standard basis of $\Bbb{R}\times \Bbb{R}^{n-1}\cong\Bbb{R}^n$) is
\begin{align}
[D\zeta_{(t,\theta)}]&=
\begin{pmatrix}
\xi(\theta) & t\frac{\partial \xi}{\partial \theta^1}(\theta) & \cdots &
t\frac{\partial \xi}{\partial \theta^{n-1}}(\theta)
\end{pmatrix}
\end{align}
So, the gram matrix of $\zeta$ looks like
\begin{align}
[\zeta'_{ij}]&=
\begin{pmatrix}
\left\langle \xi(\theta),\xi(\theta)\right\rangle & \cdots&
\left\langle \xi(\theta),t\frac{\partial \xi}{\partial \theta^i}(\theta)\right\rangle & \cdots\\
\vdots &\qquad & \\
\left\langle \xi(\theta),t\frac{\partial \xi}{\partial \theta^i}(\theta)\right\rangle & \qquad &\bigg[\left\langle t\frac{\partial \xi}{\partial \theta^{a}}(\theta),t\frac{\partial \xi}{\partial \theta^{b}}(\theta)\right\rangle\bigg]_{a,b\in\{1,\dots, n-1\}}& \qquad\\
\vdots
\end{pmatrix}\\\\
&=\begin{pmatrix}
1 & O\\
O & \qquad [(t\xi)'_{a,b}]_{a,b\in\{1,\dots, n-1\}}
\end{pmatrix}
\end{align}
where $O$ are zeros of appropriate size. This is because $\xi(\theta)$ lies on $S^{n-1}$, hence has unit norm, thus the $1$ in the top left corner; we have all those zeros, because they're inner products of $\xi(\theta)$ with the vectors $t\frac{\partial \xi}{\partial \theta^i}(\theta)$, i.e its the inner product of a vector normal to the sphere with something tangent to the sphere. Therefore, the determinant of the gram matrix for $\zeta$ reduces to the determinant of the gram matrix for $t\xi$, i.e the parametrization of the sphere of radius $t$. Hence,
\begin{align}
g(r)&=\int_0^r\int_Uf(x+t\xi(\theta))\cdot\sqrt{\det \zeta'_{ij}}\,d\theta\,dt\\
&=\int_0^r\int_Uf(x+t\xi(\theta))\cdot\sqrt{\det (t\xi)'_{ab}}\,d\theta\,dt\\
&=\int_0^r\int_{\partial B_t(x)}f(y)\,d\sigma_t(y)\,dt
\end{align}
This completes the proof.
Note that in this case we were aided heavily by the orthogonality; otherwise there would be some extra factors. More generally, we have the following:
Let $(X,g_1),(Y,g_2),(Z,g_3)$ be three Riemannian manifolds, $\phi:X\times Y \to Z$ a smooth map which is a diffeomorphism onto its image, and suppose for each $x\in X$, we have $\phi_x:=\phi(x,\cdot)$ is a diffeomorphism of $Y$ onto its image, which we denote $Y_x$ (an embedded submanifold of $Z$). Let $\lambda_{g_1},\lambda_{g_2},\lambda_{g_3}$ denote the induced Riemann-Lebesgue measures on the manifolds $X,Y,Z$ respectively. Also, for each $x\in X$, let $\iota_x:Y_x\to Z$ be the canonical inclusion map. Then, for any Borel-measurable $f:Z\to [0,\infty]$, and any Borel-measurable $E\subset X\times Y$, we have
\begin{align}
&\quad\int_{\phi(E)}f(z)\,d\lambda_{g_3}(z)\\
&=\int_X\int_{\phi_x[E_x]}f(z)\frac{d\lambda_{\phi^*(g_3)}}{d(\lambda_{g_1}\otimes \lambda_{g_2})}(x,\phi_x^{-1}(z))\cdot \frac{d\lambda_{(\phi_x^{-1})^*(g_2)}}{d\lambda_{\iota_x^*(g_3)}}(z)\,d\lambda_{\iota_x^*(g_3)}(z)\,d\lambda_{g_1}(x),
\end{align}
where the derivatives are the Radon-Nikodym derivatives of the various measures involved, and $E_x:=\{y\in Y\,:\, (x,y)\in E\}$ denotes the $x$-cross-section of $E$ in $Y$.
What this is saying is that if we have a family of parametrized submanifolds $Y_x:= \phi_x[Y]$ inside a given manifold $Z$, then we can integrate over each submanifold, and then integrate over $X$ (the space of parameters). In the above case, we had $X$ being an open interval and $Y$ being the unit sphere and $Z$ being $\Bbb{R}^n$.
To make the notation of the theorem more manageable, let $\xi,\eta,\zeta$ denote the measures $\lambda_{g_1},\lambda_{g_2},\lambda_{g_3}$ respectively, and let $\sigma_x:=\lambda_{\iota_x^*(g_3)}$ be the surface measure induced on the submanifold $Y_x$. Also, let $J_{\phi}$ denote the Radon-Nikodym derivative $\frac{d\lambda_{\phi^*(g_3)}}{d(\lambda_{g_1}\otimes \lambda_{g_2})}$; this is a smooth function $X\times Y\to \Bbb{R}^+$, and lastly, let $J_{\phi_x^{-1}}$ denote the Radon-Nikodym derivative $\frac{d\lambda_{(\phi_x^{-1})^*(g_3)}}{d\lambda_{\iota_x^*(g_3)}}$; this is a smooth function $Y_x\to \Bbb{R}^+$.
These factors are like the Jacobian determinants (with respect to the appropriate metrics) of the full mapping $\phi$ and the partial mapping $\phi_x^{-1}$. Now, the theorem can be written as
\begin{align}
\int_{\phi[E]}f(z)\,d\zeta(z)&=\int_X\int_{\phi_x[E_x]}f(z)\cdot J_{\phi}(x,\phi_x^{-1}(z))\cdot J_{\phi_x^{-1}}(z)\,d\sigma_x(z)\,d\xi(x).
\end{align}
So, this is saying that to integrate, we have to take into account the full Jacobian as a result of $\phi$ being a mapping $X\times Y\to Z$, and then if we want to pass to the family of submanifolds $Y_x$, we need the extra factor of the Jacobian determinant arising from the partial mapping $\phi_x^{-1}:Y_x\to Y$
One final remark: in the case of the ball being parametrized by the interval and unit sphere, you can verify that these Jacobian determinant factors are $t^{n-1}$ and $\frac{1}{t^{n-1}}$ respectively, and hence their product simplifies nicely, which is why we didn't have any other factors in the above example. Geometrically, this happened because of the orthogonality of the interval and the unit sphere, as I've briefly alluded to above.
