How to define smooth random initial conditions consistent with boundary conditions for this 1+2D PDE?

I wrote a fully working code to numerically solve a non-linear variation of the classic 1+2D waves equation. Here’s a stripped down version of it (classic waves equation), with static initial conditions (initial velocities are all 0):

Clear["Global`*"]

size = 30;
tmax = 20;

equDiff = D[phi[t,x,y],t,t] - D[phi[t,x,y],x,x] - D[phi[t,x,y],y,y]) == 0;

bones = Table[{x,y,RandomReal[{-1,1}]}, {x,-size,size,2}, {y,-size,size,2}];
flesh = Interpolation[Flatten[bones, 1], Method -> "Spline"];

Phi0[t_, x_, y_] := flesh[x,y] (* Random function of X and Y only. *)

boundariesCond = {
   phi[t,-size,y] == Phi0[t,-size,y],
   phi[t,size,y] == Phi0[t,size,y],
   phi[t,x,-size] == Phi0[t,x,-size],
   phi[t,x,size] == Phi0[t,x,size]
   };

initialCond = {phi[0,x,y] == Phi0[0,x,y], (D[phi[t,x,y],t]/.t -> 0) == D[Phi0[t,x,y],t]/.t -> 0};

fieldSolution = NDSolve[
   Flatten@{equDiff, initialCond, boundariesCond}, phi,
   {t,0,tmax}, {x,-size,size}, {y,-size,size},
   Method -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid","MaxPoints" -> 50,"MinPoints" -> 50,"DifferenceOrder" -> 2}}
   ];

Manipulate[
 Plot3D[
  Evaluate[phi[t,x,y]/.fieldSolution/.t -> time],
  {x,-10,10}, {y,-10,10},
  PlotPoints -> ControlActive[20, 60],
  PlotRange -> {{-10,10}, {-10,10}, {-3,3}},
  Axes -> True,
  Mesh -> True,
  MeshFunctions -> {(#3&)},
  MeshStyle -> GrayLevel[0.25],
  ColorFunction -> "Rainbow",
  ImageSize -> 500,
  Method -> {"RotationControl" -> "Globe"},
  SphericalRegion -> True
  ],
 {time, 0, tmax, 0.1,ImageSize -> Large,Appearance -> {"Labeled", "Closed"},AppearanceElements -> {"InputField", "Slider"}},
 ControlPlacement -> Bottom
 ]

This code is working pretty well on my system.

My problem is to define smooth (and derivable) random initial conditions for both the initial field and its initial time derivative (or "velocities"), that would be consistent with the boundary conditions. Currently, the code above is using only static initial conditions, so all initial velocities are 0. This isn’t very "natural" to me. So how would you define smooth and derivable random initial conditions for this system, that would be compatible with this kind of boundary conditions?

All the code should remain the same (unless you notice some mistakes or problems), except the bones and flesh functions that would need to change, so Phi0[t_, x_, y_] := flesh[x,y] would be replaced by Phi0[t_, x_, y_] := someNewRandomFunction[t,x,y].