We describe the parallel implementation of the FEMAO model for an ice-covered sea with 2D Hilbert-curve domain decomposition. Load balancing is crucial because performance depends on the local depth. We propose, compare, and discuss four approaches to load balancing. The parallel library allowed us to modify the original sequential algorithm as little as possible. The performance increases almost linearly (tested with up to 996 CPU cores) for the model of the shallow White Sea.