Abstract. Sonar surveys provide an effective mechanism for mapping seabed methane flux emissions, with Arctic submerged permafrost seepage having great potential to significantly affect climate. We created in situ engineered bubble plumes from 40 m depth with fluxes spanning 0.019 to 1.1 L s −1 to derive the in situ calibration curve (Q(σ )). These nonlinear curves related flux (Q) to sonar return (σ ) for a multibeam echosounder (MBES) and a single-beam echosounder (SBES) for a range of depths. The analysis demonstrated significant multiple bubble acoustic scattering -precluding the use of a theoretical approach to derive Q(σ ) from the product of the bubble σ (r) and the bubble size distribution where r is bubble radius. The bubble plume σ occurrence probability distribution function ( (σ )) with respect to Q found (σ ) for weak σ well described by a power law that likely correlated with small-bubble dispersion and was strongly depth dependent. (σ ) for strong σ was largely depth independent, consistent with bubble plume behavior where large bubbles in a plume remain in a focused core.(σ ) was bimodal for all but the weakest plumes. Q(σ ) was applied to sonar observations of natural arctic Laptev Sea seepage after accounting for volumetric change with numerical bubble plume simulations. Simulations addressed different depths and gases between calibration and seep plumes. Total mass fluxes (Q m ) were 5.56, 42.73, and 4.88 mmol s −1 for MBES data with good to reasonable agreement (4-37 %) between the SBES and MBES systems. The seepage flux occurrence probability distribution function ( (Q)) was bimodal, with weak (Q) in each seep area well described by a power law, suggesting primarily minor bubble plumes. The seepage-mapped spatial patterns suggested subsurface geologic control attributing methane fluxes to the current state of subsea permafrost.