© 2020 Elsevier B.V.The proposed method aims to approximate a solution of a fluid–fluid interaction problem in case of low viscosities. The nonlinear interface condition on the joint boundary allows for this problem to be viewed as a simplified version of the atmosphere–ocean coupling. Thus, the proposed method should be viewed as potentially applicable to air–sea coupled flows in turbulent regime. The method consists of two key ingredients. The geometric averaging approach is used for efficient and stable decoupling of the problem, which would allow for the usage of preexisting codes for the air and sea domain separately, as “black boxes”. This is combined with the variational multiscale stabilization technique for treating flows at high Reynolds numbers. We prove the stability and accuracy of the method, and provide several numerical tests to assess both the quantitative and qualitative features of the computed solution.