We can try to make a model with some basic (narrowly applicable) assumptions:
- Bubble is Ideal Gas
- Isothermal expansion of the bubble as it rises
- The bubble remains spherical
- Quadratic Drag
- Pressure inside the bubble is equivalent to outside the bubble
- Neglecting pressure gradient across bubble as it expands
Newtons Second Law:
$$ F_b - W - F_D = m \ddot h \tag{Eq 1}$$
For the buoyant force:
$$ F_b = \rho g V$$
Using the Ideal Gas Law (assumption 1) (isothermal expansion) to find the volume of the bubble:
$$ P_{abs} = \frac{mRT}{V} $$
Equate that to the absolute hydrostatic pressure in the fluid at depth ##h_o -h## (assumption 5):
$$ \frac{mRT}{V} = \rho g (h_o - h) + P_{atm} \implies V = \frac{mRT}{\rho g (h_o - h) + P_{atm} }$$
Substitute:
$$ F_b = \rho g V = \rho g \frac{mRT}{\rho g (h_o - h) + P_{atm} }$$
Force of weight:
$$W = mg$$
Force of Drag:
$$F_D = \frac{1}{2} C \rho \pi r^2 \dot h ^2$$
The cross-sectional area of the bubble will ultimately be a function of the depth, through the volume of the bubble ##V## and spherical bubble (assumption 3):
$$\frac{4}{3}\pi r^3 = V = \frac{mRT}{\rho g (h_o - h) + P_{atm} }$$
This implies that:
$$ r^2 = \left( \frac{3mRT}{4 \pi} \right)^{2/3} \left( \rho g (h_o - h) + P_{atm} \right)^{-2/3}$$
Now just make the following substitutions:
$$ \lambda = \rho g (h_o - h) + P_{atm}$$
$$ \implies \dot h^2 = \left( \frac{1}{\rho g} \right)^2 \dot \lambda^2$$
$$ \implies \ddot h =\frac{-1}{\rho g} \ddot \lambda$$
Making all the substitutions into ##( \rm{Eq 1})##:
$$ \frac{m}{\rho g} \ddot \lambda + \rho g m R T \lambda^{-1} -mg - \frac{1}{2} C \frac{\pi}{\rho g^2} \left( \frac{3mRT}{4 \pi} \right)^{2/3} \lambda^{-2/3} \dot \lambda^2 = 0$$That’s at least half the battle. Then just decide on some parameters, initial conditions and solve (numerically) the non-linear, second order ODE.