Save in the columns of two matrices $B$ and $C$ all vectors $b_j$ to which you applied the matrix in the previous iterations and the results $c_j=Ab_j$.
For each new system $(A+D)x'=b'$ (or $Ax=b'$, which is the special case $D=0$), approximately solve the overdetermined linear system $(C+DB)y\approx b'$, e.g., by selecting a subset of the rows (possibly all) and using a dense least square method. Note that only the selected part of $C+DB$ needs to be assembled; so this is a fast operation!
Put $x_0=By$. This is a good initial approximation with which to start the iteration for solving $(A+D)x'=b'$. In case further systems must be processed, use the matrix vector products in this new iteration to extend the matrices $B$ and $C$ on the resulting subsystem.
If the matrices $B$ and $C$ do not fit into main memory, store $B$ on disk, and select the subset of rows in advance. This allows you to keep in core the relevant part of $B$ and $C$ needed to form the least squares system, and the next $x_0$ can be computed by one pass through $B$ with little use of core memory.
The rows should be selected in such a way that they approximately correspond to a coarse discretization of the full problem. Taking five times more rows than the total number of expected matrix vector multiplies should be enough.
Edit: Why does this work? By construction, the matrices $B$ and $C$ are related by $C=AB$. If the subspace spanned by the columns of $B$ contains the exact solution vector $x'$ (a rare but simple situation) then $x'$ has the form $x'=By$ for some $y$. Substituting this into the equation defining $x'$ gives the equation $(C+DB)y= b'$. Thus in this case, the above process gives as starting point $x_0=By=x'$, which is the exact solution.
In general, one cannot expect $x'$ to lie in the column space of $B$, but the starting point generated will be the point in this cloumn space closest to $x'$, in a metric determined by the selected rows. Thus it is likely to be a sensible approximation. As more systems are processed, the column space grows and the approximation will be likely to improve a lot, so that one can hope to converge in fewer and fewer iterations.
Edit2: About the subspace generated: If one solves each system with a Krylov method, the vectors used to get the starting point for the second system span the Krylov subspace of the first right hand side. Thus one gets a good approximation whenever this Krylov subspace contains a vector close to the solution of your second system.
In general, the vectors used to get the starting point for the $(k+1)$st system span a space containing the Krylov subspace of the first $k$ right hand sides.