You can just solve Poisson's law or Gauss'law (its integral version) depending on which is the easiest as you have the complete description of the field alredy.
I think the best starting point is actually Poisson's:
$$\nabla\cdot\vec{E}={\rho\over\epsilon_0}$$
then given a field $$\vec{E}(x, y, z)$$ you get
$$\partial_x E_x+\partial_yE_y+\partial_zE_z = {\rho(x, y, z)\over\epsilon_0}$$
Whether that is easy to solve or not depends on the details.
If your surface is spherical it might be convenient to write everything is spherical coordinates.
Let's do an easy example:
$$E(r)={1\over 4\pi \epsilon_0 }{Q\over r^2} \hat{r}$$
if $r>R$ and $0$ otherwise, which is the field of a uniform spherical charge distribution over a sphere of radius $R$.
Then the divergence of the field is (in spherical coordinates)
$$\nabla\cdot\vec{E} = {Q\over 4\pi\epsilon_0} \nabla\cdot {\hat{r}\over r^2}$$
if $r>R$ and 0 otherwise. The other terms of the divergence (in $\theta$ and $\phi$, the polar angles) vanish as they are constant - the field only varies along $\hat{r}$.
So we get, using Gauss' law
$$\rho(r)=0$$ if $r<R$ and otherwise
$$\rho(r)=\epsilon_0 \nabla\cdot\vec{E}= {Q\over 4\pi} \nabla\cdot {\hat{r}\over r^2}$$
so we only need to compute
$$\nabla\cdot {\hat{r}\over r^2}$$
in $[R, \infty]$.
Now this requires a bit of calculus which I have to hand-wave [I am happy if somebody solves it more generally]
In general, $$\nabla\cdot {\hat{r}\over r^2} = 4\pi\delta(r)$$ [you can see why here] for example. However in our case we don't have a discontinuity not in the origin but in $r=R$ and I am pretty sure the result in our case would be that the divergence is $0$ everywhere and infinite in $r=R$ so
$$\nabla\cdot {\hat{r}\over r^2} = \delta(r-R)/r^2$$
meaning that
$$\rho(r)= {Q\over 4\pi r^2}\delta(r-R)$$
which is a spherical distribution of total charge $Q$ on a shell centered at the origin of radius $R$.
The main problem to deal with is the behavior of the divergence on the surface, where the charge distribution is discontinous. But in principle
$$\rho(r)=\epsilon_0 \nabla\cdot\vec{E}$$
is the solution you are looking for if you already know the field $\vec{E}$ everywhere. It's just a matter of doing derivatives.