The quantity you are dealing with is the local time of the Brownian motion at one point. It is rather difficult to deal with local times usely.
What you want to compute is the distribution of the quantity
$$\lambda(\vec r_0)=\int_0^t\delta^{(d)}\left(\vec r_0-\vec B_u\right)\mathrm du.\tag{1}$$
Note that $\lambda$ is a spatial density of time, it has dimension of $[L^{-d}T]$, where $d$ is the dimension of space.
First compute the characteristic function of the Brownian motion starting at $\vec r=0$, at time $t>0$
$$\phi(\vec q,t)=\left\langle\exp\left[\mathrm i\vec q\cdot\vec B_t\right]\right\rangle=
\int\mathrm e^{\mathrm i\vec q\cdot \vec r}\frac{\mathrm e^{-\vec r^2/4Dt}}{\left(4\pi Dt\right)^{d/2}}\mathrm d\vec r=\frac{\mathrm e^{-Dt\vec q^2}}{\left(4\pi Dt\right)^{d/2}}.\tag{2}$$
Now, we transform the delta-function in (1) into a Fourier integral thanks to the equation (2) and compute the average of $\lambda(\vec r_0)$
$$\left\langle\lambda(\vec r_0)\right\rangle=\left\langle\int_0^t\int\frac{\mathrm d\vec q}{(2\pi)^d}\mathrm e^{-\mathrm i\vec q\cdot(\vec r_0-\vec B_u)}\mathrm du\right\rangle
=\int_0^t\int\frac{\mathrm d\vec q}{(2\pi)^d}\mathrm e^{-\mathrm i\vec q\cdot \vec r_0}\phi(\vec q,u)\mathrm du.$$
The average is thus
$$\left\langle\lambda(\vec r_0)\right\rangle
=\int_0^t\frac{\mathrm e^{-\vec r_0^2/4Du}}{\left(4\pi Du\right)^{d/2}}\mathrm du.$$
There are exact results in one two and three dimensions for the average.
In one dimension, we have
$$\left\langle\lambda(\vec r_0)\right\rangle=\sqrt{\frac t{\pi D}}\mathrm e^{-\vec r_0^2/4Dt}-\frac{|\vec r_0|}{2D}\mathrm{erfc}\left(\frac{|\vec r_0|}{2\sqrt{Dt}}\right);$$
in two dimensions
$$\left\langle\lambda(\vec r_0)\right\rangle=\frac1{4\pi D}\Gamma\left(0,\frac{\vec r_0^2}{4Dt}\right)$$
($\Gamma$ is the incomplete gamma function)
and in three dimensions
$$\left\langle\lambda(\vec r_0)\right\rangle=\frac{1}{4\pi D|\vec r_0|}\mathrm{erfc}\left(\frac{|\vec r_0|}{2\sqrt{Dt}}\right)$$
I have used the same method in a recent paper
for a Brownian bridge.
To compute $\left\langle\lambda(\vec r_0)^2\right\rangle$ use the same computation technique (also described in the paper). You need to evaluate
$$\int_0^t\mathrm du\int_0^{t-u}\mathrm du'\int\frac{\mathrm d\vec q}{(2\pi)^d}\frac{\mathrm d\vec q'}{(2\pi)^d} \mathrm e^{-\mathrm i\vec q\cdot \vec r_0}\phi(\vec q,u)\phi(\vec q',u').$$
But the distribution you are looking for is not defined in dimensions larger than one: a regularization length is needed. This comes from the properties of the Brownian motion.
In the case of the Browian (the case studied in the paper), once it is regularized, the distribution is as you said: there is a $\delta$-function plus a decreasing function with the distance. In three dimensions, the decreasing function is a simple exponential.