I don't think this is possible for general $\epsilon$, and I doubt it's possible for remainder $0.0001$.
Below are some solutions with remainder less than $0.01$. I produced them by randomized search from two different initial configurations. In the first one, I only placed the circle with curvature $2$ in the centre and tried placing the remaining circles randomly, beginning with curvature $12$; in the second one, I prepositioned pairs of circles that fit in the corners and did a deterministic search for the rest.
The data structure I used was a list of interstices, each in turn consisting of a list of circles forming the interstice (where the lines forming the boundary of the square are treated as circles with zero curvature). I went through the circles in order of curvature and for each circle tried placing it snugly in each of the cusps where two circles touch in random order. If a circle didn't fit anywhere, I discarded it; if that decreased the remaining area below what was needed to get up to the target value (in this case $0.99$), I backtracked to the last decision.
I also did this without using the circle with curvature $2$. For that case I did a complete search and found no configurations with remainder less than $0.01$. Thus, if there is a better solution in that case, it must involve placing the circles in a different order. (We can always transform any solution to one where each circle is placed snugly in a cusp formed by two other circles, so testing only such positions is not a restriction; however, circles with lower curvature might sit in the cusps of circles with higher curvature, and I wouldn't have found such solutions.)
For the case including the circle with curvature $2$, the search wasn't complete (I don't think it can be done completely in this manner, without introducing further ideas), so I can't exclude that there's are significantly better configurations (even ones with in-order placement), but I'll try to describe how I came to doubt that there's much room for improvement beyond $0.01$, and particularly that this can be done for arbitrary $\epsilon$.
The reasons are both theoretical and numerical. Numerically, I found that this seems to be a typical combinatorial optimization problem: There are many local minima, and the best ones are quite close to each other. It's easy to get to $0.02$; it's relatively easy to get to $0.011$; it takes quite a bit more optimization to get to $0.01$; and beyond that practically all the solutions I found were within $0.0002$ or so of $0.01$. So a solution with $0.0001$ would have to be of a completely different kind from everything that I found.
Now of course a priori there might be some systematic solution that's hard to find by this sort of search but can be proved to exist. That might conceivably be the case for $0.0001$, but I'm pretty sure it's not the case for general $\epsilon$. To prove that it's possible to leave a remainder less than $\epsilon$ for any $\epsilon\gt0$, one might try to argue that after some initial phase it will always be possible to fit the remaining circles into the remaining space. The problem is that such an argument can't work, because we're trying to fill the rational area $1$ by discarding rational multiples of $\pi$ from the total area $\pi^3/6$, so we can't do it by discarding a finite number of circles, since $\pi$ is transcendental.
Thus we can never reach a stage where we could prove that the remaining circles will exactly fit, and hence every proof that proves we can beat an arbitrary $\epsilon$ would have to somehow show that the remaining circles can be divided into two infinite subsets, with one of them exactly fitting into the remaining gaps. Of course this, too, is possible in principle, but it seems rather unlikely; the problem strikes me as a typical messy combinatorial optimization problem with little regularity.
A related reason not to expect a clean solution is that in an Apollonian gasket with integer curvatures, some integers typically occur more than once. For instance, one might try to make use of the fact that the curvatures $0$, $2$, $18$ and $32$ form a quadruple that would allow us to fill an entire half-corner with a gasket of circles of integer curvature; however, in that gasket, many curvatures, for instance $98$, occur more than once, so we'd have to make exceptions for those since we're not allowed to reuse those circles. Also, if you look at the gaskets produced by $0$, $2$ and the numbers from $12$ to $23$ (which are the candidates to be placed in the corners), you'll find that the fourth number increases more rapidly than the third; that is, $0$, $2$ and $18$ lead to $32$, whereas $0$ $2$ and $19$ already lead to $(\sqrt2+\sqrt{19})^2\approx33.3$; so not only can you not place all the numbers from $12$ to $23$ into the corners (since only two of them fit together and there are only four corners), but then if you start over with $24$ (which is the next number in the gasket started by $12$), you can't even continue with the same progression, since the spacing has increased. The difference would have to be compensated by the remaining space in the corners that's not part of the gaskets with the big $2$-circle, but that's too small to pick up the slack, which makes it hard to avoid dropping several of the circles in the medium range around the thirties.
My impression from the optimization process is that we're forced to discard too much area quite early on; that is, we can't wait for some initial irregularities to settle down into some regular pattern that we can exploit. For instance, the first solution below uses all curvatures except for the following: 3 4 5 6 7 8 9 10 11 16 17 20 22 25 30 31 33 38 46 48 49 52 53 55 56 57 59 79 81 94 96 101 106 107 108 113 125 132. Already at 49 the remaining area becomes less than would be needed to fill the square. Other solutions I found differed in the details of which circles they managed to squeeze in where, but the total area always dropped below $1$ early on. Thus, it appears that it's the irregular constraints at the beginning that limit what can be achieved, and this can't be made up for by some nifty scheme extending to infinity. It might even be possible to prove by an exhaustive search that some initial set of circles can't be placed without discarding too much area. To be rigorous, this would have to take a lot more possibilities into account than my search did (since the circles could be placed in any order), but I don't see why allowing the bigger circles to be placed later on should make such a huge difference, since there's precious little wiggle room for their placement to begin with if we want to fit in most of the ones between $12$ and $23$.
So here are the solutions I found with remainder less than $0.01$. The configurations shown are both filled up to an area $\gtrsim0.99$ and have a tail of tiny circles left worth about another $0.0002$. For the first one, I checked with integer arithmetic that none of the circles overlap. (In fact I placed the circles with integer arithmetic, using floating-point arithmetic to find an approximation of the position and a single iteration of Newton's method in integer arithmetic to correct it.)
The first configuration has $10783$ circles and was found using repeated randomized search starting with only the circle of curvature $2$ placed; I think I ran something like $100$ separate trials to find this one, and something like $1$ in $50$ of them found a solution with remainder below $0.01$; each trial took a couple of seconds on a MacBook Pro.
The second configuration has $17182$ circles and was found by initially placing pairs of circles with curvatures $(12,23)$, $(13,21)$, $(14,19)$ and $(15,18)$ touching each other in the corners and tweaking their exact positions by hand; the tweaking brought a gain of something like $0.0005$, which brought the remainder down below $0.01$. The search for the remaining circles was carried out deterministically, in that I always tried first to place a circle into the cusps formed by the smaller circles and the boundary lines; this was to keep as much contiguous space as possible available in the cusps between the big circle and the boundary lines.
I also tried placing pairs of circles with curvatures $(13,21)$, $(14,19)$, $(15,18)$ and $(16,17)$ in the corners, but only got up to $0.9896$ with that.
Here are high-resolution version of the images; they're scaled down in this column, but you can open them in a new tab/window (where you might have to click on them to toggle the browser's autoscale feature) to get the full resolution.
Randomized search:
With pre-placed circles:
This answer requires some spatial examples, so I think it will be better if I generalize the case of a simple cube, which can be applied to all possible situations.
As you said, the total area covered by the paint in a cube with lenght $a$ is equal to:
$$
S_{\mathrm{ext}}=6\cdot a^2
$$
But then the division into $27$ cubes with lenght $1/3$ of the original gives the new surface area, greater than the original because new slices prodce more internal surface; so now the total surface is expressed by:
$$
S_{\mathrm{tot}}=27\cdot 6\cdot (a/3)^2=\frac{27\cdot 6}{9}\cdot a^2=3\cdot 6\cdot a^2=3S_{\mathrm{ext}}
$$
Now think about a Rubik's Cube, with only the cubes facing the outside are painted: since you know that the total amount of paint occupies the external surface $S_{\mathrm{ext}}$, then to know the surface not covered in paint you have to simply subtract to the total sliced surface $S_{\mathrm{tot}}$ the external one:
$$
S_{\mathrm{unpaint}}=S_{\mathrm{tot}}-S_{\mathrm{ext}}=3S_{\mathrm{ext}}-S_{\mathrm{ext}}=2S_{\mathrm{ext}}=12\cdot a^2
$$
Now you have basically derived the sum of the internal surface of every little cube, because the sum of the external is nothing more than the total covered in paint, or $S_{\mathrm{ext}}$ (with the Rubik's Cube case, the only cubes coloured).
EDIT: If the number of slices is not defined, it is generalized to a single parameter $n\in \mathbb{N}$, then the new surface tends to be greater as $n$ gets larger.
Let's assume that one cut is repeated on every edge of the cube, so if an edge is divided in three (with two slices), then the number of cubes created is expressed as follows:
$$
N_{\mathrm{cubes}}=(n+1)^3
$$
So if I "slice" a cube for every edge once ($n=1$), I get $N_{\mathrm{cubes}}=(2)^3=8$ little cubes.
The surface of each sub-cube is expressed by:
$$
S_{\mathrm{sub-cb}}=6\cdot\left(\frac{a}{n+1}\right)^2
$$
Then the total surface is expressed as th single one multiplied by the total number of sub-cubes created $N_{\mathrm{cubes}}$:
$$
S_{\mathrm{tot}}=(n+1)^3\cdot 6\cdot\left(\frac{a}{n+1}\right)^2=(n+1)\cdot 6\cdot a^2=(n+1)\cdot S_{\mathrm{ext}}
$$
Now, applying the same deduction as before, the surface not covered with paint is given by a simple subtraction between the total surface and the ponly painted one, which is the external one:
$$
S_{\mathrm{unpaint}}=S_{\mathrm{tot}}-S_{\mathrm{ext}}=(n+1)\cdot S_{\mathrm{ext}}-S_{\mathrm{ext}}=n\cdot S_{\mathrm{ext}}
$$
Now remember that the parameter $n$ is the number of slices, and not the sub-cubes created; as your previous example, the slices were $2$, and so this number appeared inside the final result.
Best Answer
As user vadim123 said in the comments, a full answer to this question appears on this page on Wikipedia.