First point: the authors tend to take a very "stringent" view of "worlds" in the many worlds interpretation. I take issue with this. I think "many worlds" is a misnomer. Better names for the Everettian interpretation would be "strict unitary evolution", "universal wavefunction" or, my favorite, "church of the large Hilbert space". For the "many worlds" are "soft" boundaries rather than "hard" boundaries. This is especially because a "splitting" of worlds can always, in principle be undone by the reverse of the unitary that split them.
Now to address your question. First of all, I would say there's no point in discussing this question for a system which is not closed. Open systems do not demonstrate energy conservation almost by definition. No one expects them to. So we can leave your possibilities 1 and 2 out of the question because they are situations that are under-specified to determine the global conservation of energy.
Response to question a): So I reiterate, this insistence on the reality of worlds in MWI is naive in my opinion. It is true that a wavefunction that might be able to be expressed as a single term might evolve into one that requires multiple terms, but it can evolve backwards too. Also, sometimes what looks like many terms might only be one term if you express it in a different basis. It might be true "for all practical purposes" that you can't reverse the wavefunction in a particular case. But when we are discussing the mathematics we don't care an ounce about practicality, and the Everettian interpretation is a very mathematical one.
So do we say energy leaves one "branch" and goes to another? No.. I find that language to be horribly confusing. Let's forget the language of branches and just consider the mathematical evolution of a wavefunction under the Hamiltonian.
First, we make a much stronger statement than "The average energy of the system is conserved". Rather, more strongly, the magnitude of the universal wavefunction* to have a given energy is conserved over time. Take state $|\psi\rangle$ and express it in the energy basis:
$$\tag{1}
U|\psi\rangle = e^{-i\frac{H}{\hbar}t} \sum_n c_n |n\rangle = \sum_n e^{-i \frac{E_n}{\hbar}t} c_n |n\rangle
$$
The Hamiltonian cannot mix states of different energy, so the weight of each term, $|c_n|^2$ is conserved.
I'm digressing a little. I want to look more at the example in the linked question of photon absorption. The atom-photon Hamiltonian in a closed system (such as a high quality optical cavity) is $(\hbar=1)$, in the rotating frame (the rotating frame makes the Hamiltonian time independent, alternatively we could consider a system that just has this Hamiltonian out the door):
$$
H = -\Delta a^{\dagger} a + \frac{\Omega}{2}(a^{\dagger}\sigma^+ + a\sigma^-)
$$
$\Delta = \omega_{\text{photon}} - \omega_{\text{atom}}$. In this rotating frame, if $\Omega = 0$ (no coupling), then we see that the state $|0, g\rangle$ has 0 energy, the state $|0, e\rangle$ also has 0 energy (this is because the choice of rotating frame). The states $|1, g\rangle$ and $|1, e \rangle$ each have $-\Delta$ energy.
An absorption transition from $|1, g\rangle \rightarrow |0, e\rangle$, naively, looks like it does not conserve energy by an amount $\Delta$.
However, such a transition does not happen unless $\Omega$ is non-zero. But, when $\Omega$ is non-zero we can see that $|1, g\rangle$ is no longer actually an eigenstate! The actual eigenstates are superpositions of the $|1, g\rangle$ and $|0, e\rangle$ states: $|+\rangle$ and $|-\rangle$, the so-called dressed states!
\begin{align}
|+\rangle =& \sin(\theta)|g\rangle + \cos(\theta)|e\rangle\\
|-\rangle =& \cos(\theta)|g\rangle - \sin(\theta)|e\rangle
\end{align}
Where $\theta$ is defined by $\tan(2\theta) = -\Omega/\Delta$. Inverting this transformation we have:
\begin{align}
|g\rangle =& \sin(\theta)|+\rangle + \cos(\theta)|-\rangle\\
|e\rangle =& \cos(\theta)|+\rangle - \sin(\theta)|-\rangle
\end{align}
So if the system starts in $|1, g\rangle$ it is actually already a superposition of $|+\rangle$ and $|-\rangle$, so it ALREADY contains a superposition of energy states. The energies are given by
$$
E_{\pm} = -\frac{\Delta}{2} \pm \frac{\sqrt{\Omega^2 + \Delta^2}}{2}
$$
For $\Omega \rightarrow 0$ we have $E_+\rightarrow 0$ and $E_-\rightarrow -\Delta$. Because $|1, g\rangle$ is not an eigenstate of the system it can change, under Hamiltonian evolution, into a different state, for example $|0, e\rangle$. Under our naive thoughts about energy it looks like the system has changed in energy by amount $\Delta$. But, if you decompose the system into the actual energy eigenbasis, you will see that the weightings of the dressed states, $|\pm\rangle$ remain unchanged other than their relative phase.
Not sure what more to say here.. Eq. (1) shows us that a system that starts in a superposition of energy states always remains in the same-weighted superposition of those energy states. This is the best quantum mechanics can do. This is what I mean when I say quantum mechanics conserves energy. This of course implies the weaker statement that the average energy of a system is conserved.
The second example reminds us that when systems are coupled their energies change and we must think in the "dressed" basis if we want to be very careful about tracking energy in the system. We must recall that the undressed states are not eigenstates of the system. This example is very relevant when thinking about absorption of off-resonant photons which is what spurred the discussion in the first place.
*Commonly called the probability of that component, but of course in MWI probability is a little controversial.
