BJT: why does Ic increase when temperature increases?
I've discussed this issue at least four times here. Maybe more. It's due to the temperature variation in the saturation current for the device, which has an opposite sign and greater magnitude than the Shockley equation provides.
Shockley Equation
As I'm sure you can find repeated, in many places, the simple Shockley equation that operates reasonably well for a BJT in active mode:
$$I_\text{C}=I_{\text{SAT}\left(T\right)}\left(e^{\left[\frac{V_\text{BE}}{\eta \:V_T}\right]}-1\right)\\ \text{where the thermal voltage is }V_T=\frac{k\: T}{q}\label{eq1}\tag{Shockley equation}$$
(The emission coefficient, \$\eta\$, is usually taken to be 1, by default.)
The saturation current \$I_{\text{SAT}\left(T\right)}\$ is effectively a \$y\$-axis intercept point found by projecting the log-chart plot of collector (or base) current vs \$V_\text{BE}\$ voltage. See this chart for a little clarity about that:

(Note that Region II is the usual &quot;operation&quot; area when in active mode.)
It's also pretty obvious that the additional -1 term in the equation removes this offset and brings the results into the general expectation of zero collector current when the base-emitter junction voltage is \$V_\text{BE}=0\:\text{V}\$.
Solving the above for \$V_\text{BE}\$:
$$V_\text{BE}=\eta\:V_T\:\operatorname{ln}\left(\frac{I_\text{C}} {I_{\text{SAT}\left(T\right)}}+1\right)\approx \eta\:\frac{k\: T}{q}\:\operatorname{ln}\left(\frac{I_\text{C}} {I_{\text{SAT}\left(T\right)}}\right)$$
From the \$\ref{eq1}\$, we can get:
$$\frac{\text{d}\: I_\text{C}}{\text{d}\: T}\approx -\frac{V_\text{BE}\:I_\text{C}}{\eta\:V_T\:T}$$
Because of the sign here, you'd expect that the collector current to decline with increasing absolute temperature. (This equation only quantifies the expectation you'd get from a cursory examination of the \$\ref{eq1}\$.)
So it's a good question you ask.
Saturation Current
Perhaps you noticed my use of \$I_{\text{SAT}\left(T\right)}\$, earlier? As it turns out, the temperature dependence of \$I_{\text{SAT}\left(T\right)}\$ is huge. The approximate equation looks like this:
$$I_{\text{SAT}\left(T\right)}=I_{\text{SAT}\left(T_\text{nom}\right)}\cdot\left[\left(\frac{T}{T_\text{nom}}\right)^{3}e^{\frac{E_g}{k}\cdot\left(\frac{1}{T_\text{nom}}-\frac{1}{T}\right)}\right]$$
\$E_g\$ is the effective energy gap (in eV) and \$k\$ is Boltzmann's constant (in appropriate units.) \$T_\text{nom}\$ is the temperature at which the equation was calibrated, of course, and \$I_{\text{SAT}\left(T_\text{nom}\right)}\$ is the extrapolated saturation current at that calibration temperature.
This formula heavily depends upon fundamental thermodynamics theory and the Boltzmann factor (do not confuse this with the Bolzmann constant \$k\$), which you can easily look up and is above represented by the factor: \$e^{\frac{E_g}{k}\cdot\left(\frac{1}{T_\text{nom}}-\frac{1}{T}\right)}\$. (It's based on the simple ratio of the numbers of states at different temperatures; really no more complex than fair dice used in elementary probability theory. Perhaps the best introduction to the Boltzmann factor is C. Kittel, &quot;Thermal Physics&quot;, John Wiley &amp; Sons, 1969, chapters 1-6 in particular.)
(Note that the power of 3 used in the equation above is actually a problem, because of the temperature dependence of diffusivity, \$\frac{k T}{q} \mu_T\$. And even that, itself, ignores the bandgap narrowing caused by heavy doping. In practice, the power of 3 is itself turned into a model parameter.)
Once you have included the equation exposing the temperature-dependence of the BJT saturation current into the \$\ref{eq1}\$ above, you find that the computation of \$\frac{\text{d}\: I_\text{C}}{\text{d}\: T}\$ now has a positive sign instead of a negative one. And this is what explains the effect you noted in your question.
I'll leave the derivation of newly arriving equation for \$\frac{\text{d}\: I_\text{C}}{\text{d}\: T}\$ to you, though.