Question: How to simplify units?

How do I get Maple 2023 to simplify/combine units in results?  For example,

Quantity(3.3897859*Unit(km^2)/Unit(m^2), 0.2) ;

should simplify to

Quantity(3389785.9, 0.2) ;

The full example where the problem occurs is given below.

alias(l = log10, l100 = log[100], pi[0] = Pi, r = sqrt, S_ellipsoid = ellipsoid) :
with(ScientificErrorAnalysis) :
with(Units) :

_km := Unit(km) :
_lm := Unit(lm) :
_lx := Unit(lx) :
_m := Unit(m) :
_rev := Unit(rev) :
alias(Q = Quantity) :
AddUnit(astronomical_unit, context = astronomy, default = true, conversion = 149597870700*m) :
pi := pi[0] :
S_spheroid := (a, b) -> S_ellipsoid(a, a, b) :

_AU := Unit(AU) :
a_Mars := Q(227939366., 1.)*_km : ##
a_Terra := Q(149598023., 1.)*_km : ##
alpha_Phobos := Q(9517.58, 0.01)*_km : ##
B_Phobos := Q(0.071, 0.012) :
d_x_Phobos := Q(25.90, 0.08)*_km :
d_y_Phobos := Q(22.60, 0.08)*_km :
d_z_Phobos := Q(18.32, 0.06)*_km :
E_I := _lx :
H_Mars := -Q(1.5, 0.1) : ##
L_Sol := Q(3.75E28, 0.01E28)*_lm : ##
m_E := -Q(14.18, 0.01) : ##
pi_Phobos := Q(9234.42, 0.01)*_km : ##
r_e_Mars := Q(3396.2, 0.1)*_km :
r_e_Terra := Q(6378137.0, 0.1)*_m : ##
S_sphere := r -> S_spheroid(r, r) :
theta_rev := _rev :

a_Phobos := (alpha_Phobos + pi_Phobos)/2 :
Delta_Sol := a_Terra - r_e_Terra : 
l_A := _AU :
m := E -> m_E - 5*l100(E/E_I) :
r_x_Phobos := d_x_Phobos/2 :
r_y_Phobos := d_y_Phobos/2 :
r_z_Phobos := d_z_Phobos/2 :
theta_rev2 := theta_rev/2 :

Delta_Mars := r(a_Mars^2 + Delta_Sol^2) : 
q := theta -> 2*((sin(theta)/pi) + (1 - (theta/theta_rev2))*cos(theta))/3 :
r_e_Phobos := (r_x_Phobos + r_y_Phobos)/2 :
S_Phobos := S_ellipsoid(r_x_Phobos, r_y_Phobos, r_z_Phobos) :

Delta_Phobos := Delta_Mars : 
L_Phobos := B_Phobos*L_Sol*S_Phobos/(2*S_sphere(a_Mars)) :
mu_Mars := 5*l(a_Mars*Delta_Mars/l_A^2) : 
rho_e_Mars := arcsin(r_e_Mars/Delta_Mars) :
theta_Mars := arccos((a_Mars^2 + Delta_Mars^2 - Delta_Sol^2)/(2*a_Mars*Delta_Mars)) :

E_Phobos := L_Phobos/S_sphere(Delta_Phobos) :
m_Mars := H_Mars + mu_Mars - 5*l100(q(theta_Mars)) :
rho_e_Phobos := arcsin(r_e_Phobos/Delta_Phobos) :
rho_o_Phobos := arctan(a_Phobos/Delta_Phobos) :

Deltarho_Phobos := rho_o_Phobos - rho_e_Mars - rho_e_Phobos :
m_Phobos := m(E_Phobos) :

Deltam_Phobos := m_Phobos - m_Mars :

Deltarhostar_Phobos := Deltarho_Phobos/Deltam_Phobos :

"Phobos apparent logarithmic brightness in astronomical magnitudes" = combine(m_Phobos, errors) ;
"Mars-Phobos apparent angular separation in arcseconds" = combine(convert(Deltarhostar_Phobos, units, arcsec), errors) ;

(*
https://www.nicolesharp.net/wiki/Solar_System_data_for_Maple
https://en.wikipedia.org/wiki/astronomical_unit
https://en.wikipedia.org/wiki/Star_Sol
https://en.wikipedia.org/wiki/Planet_Terra
https://en.wikipedia.org/wiki/Terran_radius
https://en.wikipedia.org/wiki/WGS84
https://en.wikipedia.org/wiki/Planet_Mars
https://en.wikipedia.org/wiki/Satellite_Phobos
https://en.wikipedia.org/wiki/Solar_System_by_size
https://en.wikipedia.org/wiki/illuminance
*)